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ABSTRACT 


This hybrid simulation demonstrates the feasibility of using 
a digital computer for the realization of a discrete filter-con- 
troller to drive a 3"/50 gun mount in response to step, ramp and 
stabilization orders. The gun position designation signals are 
subject to random environmental noise and the observation of 
position is contaminated by random measurement noise. These 
noise sources are assulfléd'to exhibit a normal distribution with 
zero mean. A filter is used to predict plant position and velocity 
from a noisy observable position, and a controller is used to 
provide the desired feedback of these states for a ripple-free 
response. The simulation results indicate that there is no de- 
tectable difference in response between constant filter gains and 
adaptive filter gains for this time-invariant plant. The results 
also demonstrate that the application of good programming tech- 
niques is paramount in minimizing the timing and scaling 


limitations imposed by hybridization. 








Section 


75 
2 
3 
4 
D. 
6 
7 
8 


9. 
10. 
1 
IZ. 
1ο, 
14. 
Boi 
ee 
Ke 
16; 
P. 
ZU. 


Appendix 
T. 
IS: 
Tue, 


We (c x 


TABLE OF CONTENTS 


Introduction 

Description of 3"/50 Gun Mount 
Gun Position Orders 

State Space Description of the Plant 
The Controller 

The Filter 

The Filter-Controller 


Page 
1] 
ll 
16 
17 
20 
22 
24 


Description of Hybrid Simulation Equipment 26 


Signal Conditioning and Number Scaling 


Arithmetic & Matrix Subroutines 
Program Optimum Controller 
Simulation Results 

Hybrid Control System Limitations 
Selection of Sampling Interval 
Selection of Noise Values 
Programs Hybrid I & II 

Adaptive Filter Gains 
Simulation Results for Hybrid II 
Advantages of Digital Control 
Conclusions 


Bibliography 


State Space Representation of a System 
Fortran 60 Programs 


Program Optimum Controller 


29 
33 
33 
34 
34 
43 
44 
48 
o9 
53 
60 
61 
63 


64 
68 
89 


Appendix Page 


IVI Program Hybrid I 98 
V. Program Hybrid II 110 
VT, Subroutines for CDC 160 Digital Programs 122 


Figure 


Nil 
2-1 


NJI 
12-1 
12-2 
1229 
12-4 
12-5 
12-6 
NEJ 
ο] 


I 
16-2 
1 6-8 
16-4 
18-1 
18-2 
NG 


LIST OF ILLUSTRATIONS 


Filter-Controller Simulation Schematic 


Conventional and Digital Control Power 
Drives 


Discrete Flow Graph of the Plant 
Train Plant 3"/50 Gun Mount 


Signal Flow Representation of Plant with 
Deterministic Inputs 


Simulation Scheme - Optimum Controller 

Simulation Results for Optimum Controller 
Simulation Results for Optimum Controller 
Simulation Results for Optimum Controller 
Simulation Results for Optimum Controller 
Simulation Results for Optimum Controller 
Simulation Results for Optimum Controller 
Simulation Results for Optimum Controller 


Curves of Stable Gain Matrix Elements 
Versus Noise Power to Signal Power Ratio 
for Sampling Interval of one Second 


Simulation Scheme - Filter Controller 
Program Time History Graphs 
Simulation Results for Hybrid II 
Effect of Computation Time Delay 
Simulation Results for Hybrid II 
Simulation Results for Hybrid II 


Simulation Results for Hybrid II Time Scaled 
by Ten 


Page 
I 
14 


19 
19 


2 
32 
35 
36 
37 
38 
39 
40 
41 


45 
41 
90 
ol 
92 
54 
SS 


96 


Figure 


18-4 


18-5 


18-6 


I-1 


πι 
Eé 
IV-l 


Simulation Results for Hybrid II Time 
Scaled by Ten 


Simulation Results for Hybrid II Time 
Scaled by Ten 


Simulation Results for Hybrid II Time 
Scaled by Ten 


Sample Flow Graph for Second Order 
Plant 


External Clock Control 
Time Sequence of One Sample 


Flowchart for Hybrid Programs 


Page 


97 


98 


99 


67 
92 
92 
Ἵν 


LIST OF TABLES 


Table Page 
2-1 Gun Mount Performance Characteris- 
tics 15 
9-1 Comparison of Equipment Ranges and 
Equivalent Values 30 
9-2 Number Scaling 30 
II-1 Definition of Terms for Program Optcon 69 
II-2 Definition of Terms for Program Filter 75 


Symbol 
ADC 
DAC 
A/D 
D/A 


(xxxx) E 


LIST OF SYMBOLS 


Analog-to-digital converter 
Digital-to-analog converter 
Analog-to-digital conversion process 
Digital-to-analog conversion process 


XXXX is a number expressed in the octal number 
System 


l. Introduction 

In the past few years a great deal of progress has been made 
in the theoretical application of state space techniques for the 
design of sampled-data control systems. The objectives of this 
paper are to apply such theory to the design of an optimal con- 
troller and Kalman filter for the control of the 3"/50 RF Twin Gun 
Mount and to construct a hybrid simulation of the system so de- 
signed, simulating the plant on the Pace TR20 analog computer 
and representing the filter-controller by a program in the CDC 
160 digital computer. 

The analytical methods outlined in References [1] [2] 

[3] and D are applied to the design of a discrete filter- 
controller vvith constant filter and controller gains to replace 
the conventional continuous feedback control scheme. These 
references assume that the system plant may be expressed as 

a set of linear, constant coefficient, differential equations. 
References [1] and [2] outline digital computer programs that 
compute the gain matrix elements for the filter, and the feed- 
back coefficients of the controller, respectively. Reference [4] 
presents a very thorough summary of the theory as applied to 
the optimum filter-controller problem and introduces the treat- 
ment of deterministic inputs as shown in the general simulation 
schematic of Figure 1-1. 

System operation with constant filter gains is compared to 
that with adaptive filter gains, where the latter are calculated 
on-line in the digital computer. 

2. Description of 3"/50 Gun Mount 
The conventional gun Conte system uses a typical synchro- 
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amplidyne drive system as shown in Figure 2-1. Position error 
is determined by a control transformer, the error voltage is then 
demodulated, amplified, compensated and applied to the ampli- 
dyne power amplifier, The train and elevation servos are 
essentially identical. The system has a coarse control synchro 
for slewing onto the target bearing quickly and a fine control 
synchro for tracking the target smoothly. The control perform- 
ance constraints are a maximum angular velocity of 30 degrees 
per second and a maximum angular acceleration of 75 degrees 
per second squared. 

To introduce the discrete filter-controller loop, the conven- 
tional control loop vvas broken at the input to the amplidyne 
power amplifier and the output shaft of the d-c drive motor as 
shown in Figure 2-1. The breaking of the system at these points 
allows either the conventional compensation loop and control or 
the discrete compensation loop and control to be used by simply 
throwing a switch. 

Reference [5] gives the following mathematical descrip- 
tion of the gun mount open-loop transfer function for train and 
elevation as determined by frequency and transient response 
techniques: 


G(s) = K(s2 +as +b 
s(s + pj) (s + P3) (s + Ρα) (8 t cs * d) 


where K is the forward path gain, and the three simple 


(2-1) 


poles are introduced by the characteristics of the d-c amplifier, 
amplidyne, and the drive motor. The complex zeroes and poles 
are introduced by the mechanical structure of the gun mount 

which causes an open loop resonant condition at about 8 cycles 


per second. 
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This transfer function can be approximated by the follovving 
second order plant vhose time constant represents the combined 
characteristics of the drive motor, and the effective friction and 
inertia at the motor shaft: 


G(s) = K 
s(s + a) 


(2-2) 
This approximate transfer function models the real system well 
enough for the purposes of this paper. 

In the real gun control system, non-linearities exist be- 
cause of saturation and because of the limits placed on velocity 
and acceleration that were designed into this control system for 
protective measures. To simplify the simulation and to be able 
to apply linear state space theory the plant is assumed to be 


linear even though this approximate transfer function of 


equation 2-2 holds only within certain bounds. 


GUN MOUNT PERFORMANCE CHARACTERISTICS 


Train Elevation 
Transfer - 
Function G(s) = 3.0 G(s) = 6.8 
s(s t 3.5) s(s t 3.5) 
Uncompensated 
ND Greater than one 0.668 
Natural 
EN 1.732 radians/sec 2.62 radians/sec 
TABLE 2-1 


Assuming unity feedback the train transfer function would 


have a closed loop response with a damping ratio slightly 


(E 


greater than one and a natural frequency of 1.732 radians/ 
second, The system is over-damped and has a very low natural 
frequency. Consequently the response would be sluggish. In 
reference [5] , personnel at the Naval Electronics Laboratory 
outline the application of a 'dominant mode' type of compensator 
which introduces a more desirable pole location having a damping 
ratio of 0.7 and a natural frequency of ten radians/second. Using 
a sampling rate of 0.1 seconds and this 'dominant mode' com- 
pensator they were able to provide a digital control to the 3"/50 
gun mount and to obtain a ripple free response to step and ramp 
inputs. 
Jj. Gun Position Orders 

The gun, however, must in addition respond to sinusoidal 
stabilization signals of roll and pitch as well as the determinis- 
tic inputs of a step anda ramp. Furthermore these signals are 
all contaminated by noise. 

For this paper the sea motions of roll and pitch are defined 


by the following equations: 


Roll = AR*sin(w,t) (3-1) 
Pitch = AP*sin(w t) (3-2) 
where 


AR is roll amplitude in degrees 
AP is pitch amplitude in degrees 
T, is the roll period in seconds, T, = 10 seconds 


To is the pitch period in seconds, T, = 5 seconds 


and the magnitude of the roll amplitude " assumed to be twice 
that of the pitch for each sea state, 

In the hybrid simulation the stabilization signals are 
sampled and subsequently represented by the output of a zero 
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order hold. It is assumed that the random noises appearing on 
the stabilization signals are Gaussian, with zero mean and 
with standard deviations that can be determined. 

In the conventional mode of operation of the gun mount, 
stabilization signals would be applied continuously. The mag- 
nitude and direction of these signals is a function of the sea 
state. They are required to keep the gun stationary with re- 
spect to a horizontal plane. When a target is designated to the 
gun, a Step input causes the gun to slew to this position at max- 
imum velocity. Upon acquiring the target, a ramp input is re- 
quired for smooth tracking of the target. 

It is the intention of this study 

1) to use a filter to estimate.the plant state variables (gun 
position and velocity), and 

2) to use a controller to provide the right combination of 
these estimated states in forming a scalar control to obtain an 
optimum response, The implementation of these ideas, as 
depicted in Figure 1-1, will be described in the succeeding 


sections. 


4. State Space Description of the Plant 


Consider the general second order plant of equation 2-2 to 
describe the train and elevation motions of the gun. This plant 
has been assumed to be a linear, time-invariant system and 
by the methods outlined in Appendix I can be described as a 
set of first order, linear differential equations given in the 
following matrix format: 


x = Fx + Du (4-1) 
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where and 


P sn 1 DR) 
0 -a k 


The state transition matrix can be found as a solution of 
the unforced state equation in the follovving manner: 

SJE 

s x(s) - x(0) = F x(s) 

x(s) = [ sI-F i x(0) 
where I is the identity matrix 
but 

[sI-F ] -l = l/s 1/s (sta) 
0 17 (sta) 

For the continuous case, 

Φι-1ρ) 5151 ferr) - |i — (a-e 4 "trto)) /a 


4-2 
The state transition matrix for the discrete case is 
®= |1 (1 -ε-81)/8 
-aT 


T,a 
where T = tto, and T is the sampling interval 
anda a is the time constant of the plant, 

If the control u is represented as the output of a zero 
order hold and therefore is constant over the sampling interval 
then the distribution matrix can be found as a solution to the 
forced state equation for the continuous case in the following 


manner: 
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t 
DE ®(t-t]) D dt 


ο 
A= ς Li 1/s 1/s(sta) 0 dt 
t 0 1/(s+a) K 


O 
A-|k |t (üa-e8€-t9)/a)at 


4 uu 


KË galtio) ar 
ta 
where K is the plant gain. 
The distribution matrix for the discrete case is 
A =|K ς ((1-ε-31)/α)αι 
0 (4-5) 
K $ caT dt T,a,K 
Then the complete discrete state equation depicted graphically 
in Figure 4-1 is given by 
x(k+1) = x(k) + Au(k) (4-6) 
By using the digital program 'Phidel' described in 
reference [2], the transition and distribution matrices for the 
train plant shown in Figure 4-2 for a sampling interval of one 
second were found to be 
$= |1.000 07277 and A = 10.475 
0.000 0.030 O bol 
5. The Controller 
It is desired to determine a set of feedback controller gains 
to control the plant in accordance with a selected performance 


index. Reference [2] defines one possible performance index 
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for a discrete sampled-data system as follows: 
(N) = Min È | x69 9x09 * Ru? -1) | (5-1) 
u(k) K=M 
where 
Q is a nxn positive definite, symmetric, constant matrix 
N is number of stages of the process 
R is a positive constant 
u is the control 
X is a nxl vector of plant states 
For the minimum time case and for the second order plant being 
considered in this paper the cost function reduces to 


2 
` mi 


Reference [2] goes on to develop the recursive algorithms for 


J(N) = x (N) + x (5-2) 


this discrete final-value control problem as follows: 


e he n 
Ss Beer 0, P5; 2I 


pa T 
A; = - A E d (5-3) 
T 
ΔΡ Δ 
& = + AA! 
1 


ËS ας 
à NC 


where 


i = 1 to N; N the number of stages of process for 


the minimum time case is equal to the order of the system 


I is the identity matrix 
ad is a Ixn vector of feedback controller gains 


and 
P; and Y are matrices involved in the intermediate steps in 


the recursive solution for the feedback gains. 
21 


From the Fortran 60 program 'Optcon', modified to minimize 
the final states, the following feedback controller gains were 
obtained for the train plant: 

At = | -1.203 -0.343 | 
A listing of program 'Optcon' is given in Appendix II. Table 
II-1 gives the definitions for the terms for this program and in 
addition suggests values of the parameters of the cost function 
equation 4-1 for the three cases of minimization of time, min- 
imization of control effort, and minimization of both time and 
control effort. 

6. The Filter 

In the controller discussion, it is assumed that all the 
states of the dynamic system are observable states. This con- 
dition is certainly not always true. In addition, the observable 
states may be measured only at the cost of some ambiguity due 
to measurement noise. The digital filter will provide a best 
estimate of all the states by weighting the past information with 
the present observable states. This weighting is performed with 
the knowledge of the environmental and measurement noise, 
both of which, it has been assumed, have independent Gaussian 
distributions with zero mean and a standard deviation that can 
be determined in some fashion. 

Ihe recursive equations 6-1, 6-2 and 6-3 for the Kalman 
filter are given without proof. References [1] and [3] describe 
proofs and D gives a detailed discussion for a general case. 
The optimum estimate of the plant states is given by the smoothing 


equation: 


x(k/k) = S(/k-1) + alk) [ew È 20/x-1) | iù) 
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where 

X(k/k-1) = ® x(k-1/k-1) + AE Iw bi 
is the predicted value of the states based upon the last best 
estimate of the states. The expected value of the environmental 
noise is Zero because it has been assumed to have a normal 
distribution vvith Zero mean. 

z(k) = y(k) t V(k) is the current value of the noisy 
observable. O 

2/k-1) = H X(k/k-1) + E is the predicted noisy 
observable based on the last value. The expected value of 
measurement noise is zero because it has been assumed to have 
a normal distribution with zero mean, 
and 

y(k) = H x(k) is the current non-noisy value of the obser- 
vable states. 

The filter design problem is based upon the proper selection 
of the gain matrix, G(k), so that the sum of the variances or 
errors associated with the estimate of individual states is 
minimized. The following two recursive relationships are re- 
quired to calculate the filter gain matrix and the conditional 
covariance matrix: 

G(k) = È (k/k-1)HT (HP(k/k-1)HT * R) i (6-2) 
and d 

P(k/K-1) = ® [Px-1/x-2) - G&-1)HP-1/x-2) ]4f + © 

(6-3) 
where 

H= (1 0) is the observability matrix 

R=E (ra | is the covariance matrix of measurement 
noise. Since there is only one observable in the present problem 
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Ris a scalar. 

Q=AE (wew? Ja? is the covariance matrix of excitation 
or signal noise. A discrete schematic of the filter is shown in 
Figure l-l. 

A stable gain matrix of 

G= | 0.786 

e 
was obtained from the Fortran 60 program 'Filter' by providing 
it with the following parameters: 
2 

oc 0.052 

R(1,1)=0.02 
and 

P(O)= {0.02 0 

0 I 

A listing of program 'Filter' is given in Appendix II. Table 
II-2 defines the symbolic terms for this program. 

7. The Filter- Controller 

The combined filter-controller simulation is shown in 
Figure 1-1. The double lines represent the transmission of 
vectors and the single lines, scalars. Deterministic inputs 
are introduced and are represented by a vector quantity, DI(k), 
with the elements consisting of the input and its successive 
derivatives, The difference vector, B(k), between the determ- 
inistic inputs and the predicted states is multiplied by the 
optimum controller gain matrix, Al, A measure of the perform - 
ance of the system can be deduced from the B vector. When 
the B vector is very near zero, the filter is estimating states 
which are approximately et to the deterministic inputs. 
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Then the plant control becomes very small and the system has 
responded to the inputs. 

The output of the controller is a scalar control u(k) which 
is constant over the sampling interval and is given by 

u&) = AT [že - DI 0) | (7-1) 
This control must be introduced so that it identically affects 
both the plant and the filter. This is done by feeding the same 
control to the filter as the plant but delayed by one sampling 
interval to account for the difference in the sample indicies. 
The controlled plant and the controlled filter equations are 
then given by 

x(k+1) = ® x(k) + Au(k) * Aw(k) (7-2) 
and 

duet. lit Vr Avalon) (7-3) 

Since the train plant of the gun is a type one system, it 
has a steady state error in response to a ramp input. This means 
that the B vector would not go to zero unless the DI vector is 
modified. This is accomplished by introducing a DA factor. 
In this case study, the DI vector is made up of two elements: 

DI(1,1) = step + ramp*t + stabilization 

DI(2,1) = ramp 
where 

step is the magnitude of the step 

ramp is the magnitude of the ramp 
and 

stabilization is the magnitude of the stabilization signal 
over one sample interval. 


Consider Figure 7-1 and a unit ramp input. 
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X (00) - 1 DI, = 
X, (00) = 0 
U(00) = 0 

but 
X, (00) = 3U(00) - 3.5X, (00) 


ο # -3.5 
hence introduce DA factor 
X, (00) = 3DA - 3.5X, (00) 
“ DA = 3.5/3 
A general expression for the DA factor for any magnitude ramp is 
given by 


DA = 3.9 ramp s 
E (7-4) 


A Fortran 60 program, 'Hybrid II', was used to simulate the 
filter-controller design and to compare the results with the hy- 
brid simulation which is described in the following sections. A’ 
listing of this program is included in Appendix II. 

8.: Description of Hybrid Simulation Equipment 

The experimental part of this thesis was carried out in the 
Hybrid Control Laboratory; United States Naval Postgraduate 
School. Simulation of the system required the following equip- 
Men 

A. The Control Data Corporation 160 digital computer is 
an electronic computer controlled by an internally stored program 
in sequential locations. Memory capacity is 4096 12-bit words 
of magnetic core storage, with a storage cycle time of 6.4 
microseconds. Instructions are executed in one to four storage 
or memory cycles with the time required for execution varying 
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FIGURE 7-3 
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from 6.4 to 25.6 microseconds or approximately 15 microseconds 
on the average. A computer word is made up of 12 binary digits 
numbered from right to left from 0 to ll. All positive numbers 
must have a zero in bit 11; all negative numbers, a one in bit 
1l. Hence positive numbers range from (0000), to (377 7)g and 
negative numbers from (7777)g to (4000)g. 

The CDC 160 is very limited in its arithmetic computation 
capability. Twelve bit addition and subtraction is accomplished 
in two or three memory cycles using ones complement notation. 
Multiplication and division can only be accomplished by success- 
ive 12 bit addition or subtraction, respectively, resulting in 
excessive programming and execution time requirements, 

B. The Control Data Corporation 168 arithmetic unit pro- 
vides the CDC 160 computer with single or double precision, 
fixed point arithmetic operations, either integer or fraction but 
not floating point, for addition, subtraction, multiplication and 
division. However, library subroutine mulop provides single 
precision multiplication in floating point format by first mul- 
tiplying the numbers and then dividing by a scale factor which 
in effect post-shifts the answer so that the binary point is in 
the correct position. Unfortunately, the execution time for 
mulop is very long; 400 microseconds for no scaling and 550 
microseconds for scaling, 

C. The Pace TR20 analog computer is a solid state computer 
having 20 operational amplifiers with a saturation voltage of 
plus or minus ten volts. By the connection of the operate-reset 
relays of the TR20 to the same relays of the ADC-DAC, remote 
control of the analog computer is possible. The run key of the 
CDC 160 starts the digital program and operates the analog 
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computer: the master clear key resets the analog computer. 

D. The ADC-DAC is a 12 bit conversion unit that functions 
in the voltage range of Zero to minus ten and has a common minus 
five volts bias source. There are 12 channels analog-to-digital 
and four channels digital-to-analog. Approximately 100 micro- 
seconds is required for each A/D number conversion and approx- 
imately 20 microseconds for each D/A number conversion. The 
input disable jack allows external timing control of the ADC 
channels. 

E. The EAI 1110 Variplotter is a solid state x-y pen 
recorder with a sensitivity of 100 microvolts per inch. It was 
used to record the simulation results, 

9. Signal Conditioning and Number Scaling 

Since the voltage ranges of the analog computer and the 
ADC-DAC, and the number range of the CDC 160 computer are 
not compatible, a certain amount of signal conditioning and 
number scaling is required. Table 9-1 shows the ranges for 
each piece of equipment and their equivalent values. An 
Ax + B transformation on the analog voltages is performed to get 
them into the range of the ADC and thence into the digital com- 
puter. A reciprocal transformation is performed for the transfer 
in the other direction. The A part of the transformation is a 
scale factor equal to 5/8 for the output analog voltages and 
equal to 8/5 for the input analog voltages. This division of 
the output analog voltages by 8 or multiplication by 273 pro- 
vides a shifting of the binary point 3 places to the right. Hence 
all voltages being transferred to the CDC 160 computer have a 
scale factor of three binary. The B part of the transformation is 
a five volt bias factor required to convert all analog voltages 
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COMPARISON OF EQUIPMENT RANGES 
AND EQUIVALENT VALUES 


ANALOG SCALED ADC OUTPUT CORRECT OBSERVED 
VOLTAGE ANALOG  -(0.625X+Bias) DI GITAL DIGITAL 
(DECIMAL) VOLTAGE VALUE VALUE 
X .625X OCTAL OCTAL 
-8 -5.000 -0.000 4000 4003 
M -4,375 -0.625 2077 4407 
-6 -3.75 -1.250 4777 5007 
-5 -3.125 -1.875 5377 5363 
= -2.500 -2.500 5777 5773 
-3 -1.875 -3.125 6377 6363 
E 21.250 -3.750 6777 7003 
nai -0.625 -4.375 70g 7403 
0 0.000 -5.000 0000/7777 0003 
1 0.625 -5.625 0400 0377 
2 10025 -6.250 1000 0777 
3 1.875 -6.875 1400 1377 
4 2.500 -7.500 2000 1777 
5 3.125 -8.125 2400 2377 
5 3.750 -8.750 3000 2777 
7 4.375 -9,375 3400 3377 
8 5.000 -10.000 3777 3773 
TABLE 9-1 
NUMBER SCALING 
DECIMAL OCTAL BINARY SCALED OCTAL 
COMPUTER WORD 
1.000 1.000 s00/100/000/000/0 0400 
-1.203 -1.151 s00/100/110/100/1 7313 
= 243 -.257  s00/0080 20094100 /] 7650 
.020 .012  s00/000/000/101/0 0005 
. 004 .002 | s00/000/000/001/0 0001 
TABLE 9-2 
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into the range of the ADC. 

By limiting the analog voltages on the TR20 to plus or minus 
eight, taking 0.625 of each variable to be outputed to the ADC, 
and adding it to a bias of plus five volts, the output of the 
summer is in the range of Zero to minus ten volts ready for in- 
put to the ADC. To go in the other direction, five volts are 
added to the DAC output and the entire quantity is multiplied by 
1.6. Then the output of the amplifier is in the same range of the 
analog computer, See Figure 11-1 showing this signal transfor- 
mation for program Optimum Controller, 

Since bit number 11 (the 12th bit of the computer word) is 
the sign bit, then the remaining 11 bits are available to repre- 
sent a number, If the decimal point of the computer word is 


scaled three places to the right, this gives the following bit 


representation: 
XXX. XXX XXX XX 
| A 
sign 2-1 
bit SE 
Dg 


This implies that numbers in the range of +7.999 with a binary 
scale factor of three can be accomodated by the computer. 
Since a scale factor of three was selected for numbers in 
the digital computer, all constants must be changed to the 
above format. Table 9-2 shows the decimal, octal, binary and 
scaled computer word representations for a few constants. Sub- 
routine mulop has a second argument called the scale factor. 
In multiplying two numbers with scale factors of three, the 
result has a scale factor of six. Hence a three bit left shift 
of the binary point must be performed by dividing by the scale 
on 
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factor. 

Since the sign bit must be included one bit of the number is 
lost and, therefore, the smallest decimal number that can be 
represented with a binary scale factor of three is .004. However, 
because of inaccuracies in the ADC- DAC conversion as is shown 
in Table 9-1, the smallest decimal number used in the hybrid 
simulation was .02. 

10. Arithmetic & Matrix Subroutines 

Subroutine mulop was the only available arithmetic sub- 
routine at the start of this study. Therefore, before any sim- 
ulation could be attempted, basic arithmetic and matrix sub- 
routines had to be written, These subroutines are described in 
detail in Appendix VI. Error halts are provided in these sub- 
routines for both summation and product overflow. However, 
no provision is made for over-riding this overflow and using the 
full positive or negative values. These subroutines are written 
for single precision or 12 bit numbers. The scaling of the num- 
bers in the computer is flexible and can be changed by modifying 
all constants entered into the computer. The scale factor in 
cell (0027), must be adjusted, as mentioned in Section 9, to 
return the binary point to the same place after multiplication or 
division. 

11. Program Optimum Controller 

To get a feel for the hybrid simulation problem and to check 
out the arithmetic and matrix subroutines it was decided to try 
a Simulation that was simple and whose results could be easily 
verified. Such a simulation problem is the optimal discrete- 


time control system described by the following equation: 
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u(x) =A" | xo - pre | (10-1) 
where the true state vector X(k) is assumed completely observ- 
able without noise contamination. This control law, using feed- 
back coefficients obtained from program''Optcon' for the minimum 
time case, will drive the state variables to zero for the regulator 
problem of initial conditions and to some specific value for the 
controller problem of deterministic inputs. The system sim- 
ulation schematic is shown in Figure 11-1. The listing and 
method of operation of program Optimum Controller are described 
in Appendix III. 

12. Simulation Results for Program Optimum Controller 

Figures 12-1 through 12-7 are the x-y recorder plots of 
the gun train plant response for various initial conditions and 
deterministic inputs that were limited to the optimal control 
region of the phase plane. In all cases optimum control is 
effectively achieved in the minimum time of two samples. The 
finite response time indicated on the control signal u(t) is due 
to the inertial lag of the EAI Variplotter. The external timing 
control was found to be the most accurate, and was used for all 
recordings. There were no observable hybrid control problems 
with this relatively simple second order plant and limited inputs. 
However, for higher order plants, hybrid limitations may be 
encountered as discussed in the next section. 

The approximation of the roll and pitch stabilization sig- 
nals by the output of a zero order hold proved to be satisfactory. 
This approximation is used in the filter-controller simulation of 
Section 18., 

13. Hybrid Control System Limitations 
The response of a system controlled by a digital computer 
34 
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may be affected by the following general limitations of hybrid- 
ization: 

A. There is a time delay due to the serial conversion of 
analog quantities to digital numbers by the ADC, The serious- 
ness of this problem is determined by the number of simultaneous 
samples required. In program Optimum Controller, where four 
A/D conversions per sampling interval are required for this 
second order plant, this difficulty was not experienced. If, 
however, the order of the plant was increased, then this effect 
might be noticeable. In such a case an external sample and 
hold device would be necessary to sample all inputs at the same 
instant of time and hold these values until the ADC completes 
the conversion on all channels. 

B. There is a time delay due to the serial nature of the 
computation of the control in a digital computer. Because of 
this step-by-step performance of arithmetic operations, a 
finite time exists between input to the digital computer of 
samples from the plant and output of the control from the com- 
puter to the plant. This time delay is a function of computer 
speed and the complexity of the computation. This problem can 
be minimized by using a high speed digital computer with faster 
arithmetic subroutines than the CDC 160 and by efficient com- 
puter programming between the sampling and the output. 

C. There is a magnitude limitation imposed by the A/D 
or D/A conversion processes, This is a physical equipment 
limitation in this case with a range of plus or minus five volts 
as discussed in Section 9, In general this magnitude limitation 
will depend upon the specitic ADC-DAC equipment, This prob- 
lem would be greatly reduced with shaft to digital encoders, 
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for example, 

D. There isa limit in the magnitude of the floating point 
constants imposed by the computer word length. As discussed in 
Section 9, the smallest number that can be represented in the 
computer with a scale factor of three is .004. A computer with 
an 18 bit word length would allow more flexibility in the com- 
putation and entering of constants. 

E. There is a sample timing control problem in that the 
computer is required to sample the inputs at regular intervals. 
This can be done either internally by subroutine delay or ex- 
ternally by a clock. The internal time delay is not as accurate 
as the clock because the actual time delay depends upon vari- 
ations in computation time of arithmetic operations that may 
vary from cycle to cycle. This fact also makes the internal time 
delay much harder to adjust. External timing control was used 
in the hybrid simulation schemes. This type of external clock 
control is much more accurate and would allow time sharing of 
the computer with other equipments. 

14. Selection of Sampling Interval 

In the selection of a sampling interval for the digital 
filter- controller a number of factors must be considered. The 
inertia of the gun is large and the system is velocity and 
acceleration limited, Too small a sampling interval would re- 
quire controls of such large magnitude in the minimum time 
case that they may cause the system to saturate and operate 
non-linearly. On the other hand, for too large a sampling 
interval the late updating of the plant may be too slow for a 
fast target. If a number of systems are to be time shared by 
the same computer, then a long sampling interval would allow 
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more systems to be controlled by the same computer. The types 
of input signals must also be considered. For instance, in 
this case the sampling interval must be small enough in order 
that the roll and pitch stabilization signals can be well rep- 
resented by the output of a Zero order hold. The limitations of 
the equipment to be used must also be considered. In this case, 
the parameters and constants had to be vvithin the number range 
of the computer. Considering all these factors, a sampling 
interval of one second was chosen, 
15 ` Selection of Noise Values 

| The choice of the values of measurement ard environmental 
noise, and the initialization of the error covariance matrix can 
be by conjecture, actual measurement of the statistics or by 
equipment limitations. In this simulation, in order to stay vvith- 
in the number range of the computer, the selection was predi- 
cated by the latter method. It is assumed that the probability 
density functions of V(k) and W(k) are rormal with zero mean. 
Both V(k) and W(k) are one dimensional random variables. Hence 
they are scalars. The noise values for the hybrid simulation 
were selected as follows: 

A. The measurement noise is given by 
R(1,1) = E [væ] = 0.02 

The standard deviation of the measurement noise is given by 


O, = 0.1414 


If a shaft to digital encoder is used to transform the shaft angle 
output to a digital number, then the ambiguity of this transfor- 
mation determines the value of R(1,1). 

B. The covariance matrix of random disturbance due to 
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environmental or excitation noise is given by 


Q=AE |w% |å = Ae A 


where 

o = W(k), is the variance of excitation noise, 
and 

AA = 0.984 0.515 


0,515 02620 
The variance of excitation noise is calculated by assuming a 
pseudo noise to signal ratio of R(1,1)/Q(1,1) » 1. 
Therefore 
Q(1,1) = .02 =o" AAt 
From this relationship the variance and the standard deviation 
of the excitation noise are calculated to be 


"E = ,052 
W 


and 


o = 1228 
W 


For this weapon, the excitation noise can be considered to be 
a random signal due to jitter. The jitter may be caused by 
noisy transmission of gun orders, sharp variations in the sta- 
bilization or designation signals, and ambiguity in bearing 
alignment between the designator and the gun. 

C. The estimation covariance matrix for the second order 
plant is a two-by-two matrix. Pa (1,1) is the initial variance 
on position and is based on the statistics of the measurement 
noise. It is therefore reasonable to assign R(1,1) as the value 
of E: (1,1). P.(2,2) is a measure of the confidence in the 
initial velocity measurement. In this simulation, the velocity 
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of the plant is unknown and must be estimated by the filter. 
Consequently, a large value compared to P,(1,1) is chosen to 
show this uncertainty. Off-dimensional terms are assigned zero 
value since the initial state estimates are uncorrelated. The 
following error covariance matrix was used in program Hybrid II 


to calculate filter gains on-line: 


D 0.02 0 
0 1 


Figure 15-1 represents a normalized plot of the gain matrix 
elements for various pseudo noise-to-signal ratios. R(1,1) is 
a measure of the expected position measurement noise power 
and Q(1,1) is a measure of the expected position signal power. 
This plot is made by considering a pseudo noise-to- signal ratio 
in which R(1,1) is made equal to multiple values of Q(1,1). For 
a noise to signal ratio of one, the stable gain matrix elements 
are always given by 

G(1,1) = 0.68 
and 

G(2,1) = 0.45 
l6. Programs Hybrid I & II 

Programs Hybrid I and II were written in machine language 
for the CDC 160 computer to simulate the discrete digital filter- 
controller part of Figure 1-1. Hybrid I is described thoroughly 
in Appendix IV, and uses constant filter and controller gains. 
Hybrid II is described in Appendix V, and uses on-line com- 
putation of filter gains. These programs solve the controlled 
plant equation 7-2 and the controlled filter equation 7-3 to 
calculate the digital control for the weapon. Figure 16-1 shows 
the simulation schematic for the hybrid control system. 
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In the writing of the Hybrid programs every effort was made 
to minimize the hybridization limitations discussed in Section 13. 
For example, in the first edition of Hybrid II, the time between 
sample and output was 250 milliseconds, However, by the 
application of good programming techniques and by the re- 
location of certain computation segments, the delay was reduced 
to 130 milliseconds. 

Figure 16-2 shows the time history for programs Optimum 
Controller, Hybrid I and II. In all three programs the A/D 
sampling takes 820 microseconds or less and therefore this 
limitation is small compared to the computation delay times. A 
computation delay of 13.6 milliseconds did not affect the res- 
ponse of the system as observed and discussed in Section 12 for 
the simulation results of program Optimum Controller. However, 
a computation delay time of 130 to 150 milliseconds for the 
Hybrid programs does affect the response of the system as shown 
in Figures 16-3, 18-1 and 18-2. In these Figures a comparison 
of the step response of the system is made between normal gun 
operation being sampled every second, and the system time 
scaled by a factor of ten being sampled every ten seconds. 
Figure 16-4 illustrates in more detail what happens when this 
compute time ὧν is significant. The time axis of this Figure 
is expanded. If Ga 0, then the first and succeeding controls 
are based on the samples of the plant at the proper instant and 
the controls are applied to the plant at that same instant pro- 
ducing the desired optimum response. However, if tc > 0, then 
the control is applied t. seconds after each sample time. In 


the mean time the plant states are advancing without the proper 
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control. When the control is finally applied it is no longer the 


optimal control for the plant at time KT + te - 


Step = 1 
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EFFECT OF COMPUTATION TIME DELAY 
FIGURE 16-4 


The computation times for the hybrid programs are about 
ten times longer than that for program Optimum Controller. The 
Hybrid programs have an increased complexity of computation but 
the predominant reason for this longer time delay is that a very 
slow multiply subroutine is called very frequently in the matrix 
calculations. Mulop takes 550 microseconds each time it is 
called. Itis felt that with newer, higher speed computers" 
and faster arithmetic subroutines, an improvement in computation 
time of ten to one 1s feasible, Therefore the results obtained by 
time scaling the system by a factor of ten are justifiable rep- 
resentation of the type of control that could be achieved with 
a high speed computer, Because of this consideration, the 
remainder of the simulations results for the Hybrid programs 
are for the time-scaled plant. 
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17. Adaptive Filter Gains 
In program Hybrid I the steady state filter gains are entered 


as constants. In program Hybrid II the filter gains are computed 
on line for each sample by calling subroutine vf gain. In this 
adaptive filter gain case, the filter gains increase from initially 
computed values, to steady state values in approximately 
seven samples, From this instant on, the adaptive filter gains 
have reached their steady state values and both Hybrid programs 
have the same response for this time-invariant plant. The sim- 
ulation results for programs Hybrid I and II indicate that there 
is no improvement in the response by using adaptive filter gains. 
In fact duplicate pen recordings of the phase plane were taken 
for step responses with one response traced on top of the other. 
Because of this exact reproduction of responses all the simul- 
ation results shown are from one program, Hybrid II. 
18. Simulation Results for Hybrid II 

Figures 18-1 through 18-6 are pen recorder graphs of the 
gun train plant responses for various deterministic inputs of 


step, ramp and stabilization signals, and for measurement and 


environmental noise. Two points must be considered in observ- 


ing the simulation results that contain either noise and/or 
Stabilization signal. The noise is reinitialized to start at the 


top of the list for each pen recording so that each individual 
recording can be compared for that particular deterministic 
input, The sinusoidal stabilization signal is not synchronized 
to start at the same position for each pen recording, and there- 
fore, the individual time axes are not correlated. 

Figures 18-1 and 18-2 illustrate again the hybridization 
problem that occurs when too long a delay between sample and 
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output of the control exists. In contrast, the time scaled step 
responses are the desired dead-beat, optimum responses. 
Figure 18-3 shows how the addition of environmental and meas- 
urement noise affect the gun response., In Figure 18-4 a pitch 
stabilization signal.is introduced with a step input. The approx- 
imation of the pitch stabilization signal by the output of a zero 
order hold is seen to be satisfactory. In Figure 18-5 the 
requirement for the DA factor to ensure zero steady state error 
of this type one plant in response to a ramp is illustrated. 
Figure 18-6 depicts a typical situation of the weapon slewing in 
response to a step to a particular bearing, and tracking a target 
in response to a ramp. Thé hump is produced in this position 
response because the DA factor is introduced three samples 
before the ramp is ihitiated. In hind sight, a more judicious 
method of introducing this DÀ factor would be to detect the 
magnitude of the deterministic velocity input and calculate 
equation 7-4 each sample. Then, if the deterministic velocity 
vector is Zero, the DA factor would also be zero. 
19. Advantages of Digital Control 

The use of a digital filter-controller has many significant 
advantages over the conventional analog system. In cases 
where more accuracy and higher sensitivity are required, 
digital compensation techniques excel over analog methods. 
In addition, virtually noise-free transmission of data in the 
form of digital numerical codes is possible. Non-linear pro- 
gramming and self-optimizing programming techniques can be 
utilized with the digital computer to provide a degree of flex- 
ibility not heretofor achievable with continuous control systems. 

High operating speeds and efficient input/output procedures 
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have made the utilization of general purpose digital computers 
in real-time sampled data systems feasible. Time sharing of 
many equipments vvith the same computer is another important 
aspect. The use of vector-matrix notation for systems that 
can be represented by linear differential equations is ideally 
suited for digital computation, and has led to development of 
programs for analysis and control of systems. 

20. Conclusions 

The folloving statements are cited as being important con- 
clusions to this hybrid simulation control problem: 

A. The results of this hybrid simulation have demonstrated 
the feasibility and simplification in design which accompany 
the use of a general purpose computer as an integral component 
of the weapon control system. 

B. The results also demonstrate that the application of a 
discrete filter-controller to provide control of a tactical weapon 
is feasible. 

C. The results demonstrated that there was no difference 
between using constant or adaptive filter gains for this time- 
invariant weapon plant. This means that the simpler com- 
putational program of Hybrid I can be used, and thereby 
minimize the computation time required to effectively control 
the gun. 

D. The application of good programming techniques to 
minimize hybridization limitations is important in providing 
optimum digital control. 

E. The time-shared control of a number of weapons sys- 


tems serially with the performance of other functions is feasible 
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for a general purpose, high speed computer having fast input, 


output capabilities and fast arithmetic subroutines . 
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APPENDIX I 
STATE SPACE REPRESENTATION OF A SYSTEM 


1-1 General Theory 
The state of a system is defined as a set of variables which 


in conjunction with the system inputs com pletely describe the 
behaviour of the system. Consider a linear, time invariant 
system described by the second order differential equation: 
Y+ay+c= f(t) (1-1) 
The state variablesx, and x? and the plant control are defined 


as follows: 


X] = y 
Ko 5 xi - y (1-2) 
u =f(t)-c 


If f(t) is known for all t greater than or equal to to and if the 
initial states are known, then the states can be determined for 
all time t. The second order linear differential equation is 
reduced to a set of first order state equations as follows: 

m ES 

Χρ SEDIS astu 
or written in matrix notation as: 

X = Fx + Du (1-3) 
where 

x is an nxl vector of the system states 

u is an rxl vector of the system inputs or controls 

F is an nxn matrix of constants 

D is an nxr matrix of constants 


here . and 
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In general the output of the system may be a linear combin- 
ation of the system states and controls, and a new set of output 
equations can be defined: 

y= Ax + Bu 
where 

y is a qxl vector of system outputs 

A is a qxn matrix of constants 

Bis a qxr matrix of constants 

It can be verified that the continuous solution to equation 
I-l is given by 

x(t) = (t-ta) xt) + < 9(-t) D u(t) dt; (1-4) 

lo 
where 

(t-t) is called the state transition matrix and is defined 
as the inverse Laplace transform of (sI-F)”I 
and 

ς (t-t) D dt} is called the state distribution matrix 

to 
for the case when the control u is the output of a zero order hold. 


The discrete solution to equation I-4 is given by 


x(k+1) = ® x(k) tA ulk) (I-5) 
where 
$- FI 
A= te Patty} D dt) 
O 


ln 15 the sampling interval 
and 
kk es (k+1)T 
The advantages of describing an nth order linear differential 
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equation by a set of first order differential equations are: 

(i) ease of solution of these first order equations by digital 
computer programming techniques 

(ii) ease of correlation between state variables, flow graphs 
and analog computer simulation allowing unified study of sampled- 
data systems 

(iii) ease of specification of initial conditions. 
I-2. Flow Graphs 

The flow graph of a linear, time invariant system can be 
obtained from the state equation I-3 or vice versa. Figure I-l 
was drawn with the following ground rules: 

(i) The system states are selected as outputs of integrators 
and the nodes selected as states have only one branch entering 
the node. 

(ii) e is connected to xi) by a unity transmission branch 


(iii) all feedback paths are from the integrator outputs of 


e 
X], X2 etc; to the node representing xj. 
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SAMPLE FLOW GRAPH 
FOR 


SECOND ORDER PLANT 
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APPENDIX II 
FORTRAN 60 PROGRAMS 

II-1. General Description 

This Appendix contains the following Fortran 60 programs 
that were used to provide constants required for the hybrid 
simulation and also provide a means of verifying the results from 
the hybrid simulation: 

(i) Optcon 

(ii) Filter 

(iii) Hybrid II 

Each program has comment cards which explain their oper- 
ation and a sample data deck is included, Table II-1 defines 
the terms for program 'Optcon' and outlines the three cases for 
selecting a desired performance index. Table II-2 defines the 


terms for program 'Filter'. 
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DEFINITION OF TERMS FOR PROGRAM OPICON 


Variable Definition Remarks 

Nstage number of stages 

N order of system 

R control effort weighting selected by the operator 
constant in the cost for desired performance 
function 

Q state weighting matrix selected by the operator 
in the cost function for desired performance 

A nxl control distribution constant for a given time- 
vector invariant plant 

È nxn state transition constant for a given time- 
matrix invariant plant 

y intermediate variable initial value Y =0 


matrix in recursive 
solution for controller 
gains 


P intermediate variable initial value P dom Q 
matrix in recursive 
solutions for controller 


gains 

AT feedback controller initial value AT = 0 
gain vector 

Case I Minimum time or mini- Q, =I, Q* =0, R=Oin 


mization of terminal states program 'Optcon' 
For second order plant the minimum time performance index is 
J(N) - x (N) 4 x (N) 


Case II Minimum effort or fuel Qo =1,Q0=0, R=lin 
program 'Optcon' 


For second order plant the minimum control effort performance 
index is 


TON) = x2 (n) +x2 (N) +min È u2(k-1) 
! 2 u(k) k=l 


TABLE II-1 
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Case III Minimum time and effort Qus lues #0, R=1 
in program 'Optcon' 


For second order plant the minimum time and control effort 
performance index is 


JON) = min = x! (k)x(k) + "m, 
zl 


u(k) | 


Note: Q and R may be selected to satisfy different cost func- 
“Paj . Another method for selecting Qand R is given in reference 
4|. 


TABLE II-1 (cont'd) 
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-- - — ———————€——— aa a —-—  —— S — S—— 5 s——- ss —- —  - — — e im — n n n 


m saksi 

MENS Μπ sË 2 jase 2») 2) sDEL( 12) AKU è 
16M(12912)sQ(12912)sX(90C) s ITITLE(12)sFM(16)sEM(12)sHM(12312)> 
1 Y1(909)»Y2(909)»Y3(900) 


TC THIS PROGRAM UTITIZES A COST FUNCTION, JINY=MINIMUMTSUM XTINITEO*XTNY#T 
C SUM R*U(N-1)**2). AN UNLIMITED NUMBER OF ITERATIONS MAY BE MADE AT 
C A COMPUTATION RATE OF 2000 PER MINUTE AFTER THE PROGRAM HAS BEEN 

———C-——COMPItED. THE OUTPUT OF THIS” PROGRAM I'S THE FEEDBACK GAIN MATRIX —— - 
C A TRANSPOSE. THE FOLLOWING RECURSIVE EQUATIONS WERE DERIVED USING 
C DYNAMIC PROGRAMMING» 

ΜΝ. ATUKTE-UDELCTXPTK-IT*PHITZ(DELT*PTK-IT*DELTR] TTT 
C PSI(K)=PHI+DEL¥AT(K) (2) 
C PI(K)-PSIT(K) *PCK- 1)*PSI (K) -QFRXA(K)*AT(K) (3) 

— C —— —P(0T20,ATU0)80,5-PSTCOYT-0 — | 
C EQUATIONS 1» 2» AND 3 CONSTITUTE THIS PROGRAM 





——————————————— t 





E C — CALL IN DATA AND INTITTACIZE A 
DO 1111 151” 
READ 309KASE 
eee Si FORMAT (TI 00m 
READ 1 sNsNSTAGE 
1 FORMAT (8110) 
te PA TT 
2 FORMATI (4F2060)) 
READ 2 ((Q(I9J) sJ=1lsN) sI=l oN) 
E cC — PRINT, THe DATA READ ING A0,.--—-—.. OOOO 
PRINT 100 
100 FORMAT(1H1) 
Me RR SAS EH 
32 FORMAT(9Xs5HKASE= 913) 
PRINT 3sNsNSTAGE 
3-FORMAT(U/79X,2HN7,13,»20X, THNSTAGE*2,I131^ — 
PRINT 4sRsDT 
4&4 FORMAT(/9Xs2HR=sF1501192Xs3HDT=sF15% 11) 
PRENT os No (JI TI NIKI SIN om 
9 FORMAT(/9X»7HQ(I5J)2/(2F15.11) D 
CALL PHIDELI(PHI sDEL sN, DT) 
DO S EPN -’''''-------------------- S eee Gc  ÁÓ 
DO 5 Jz1»N : 
GM(I,»J)20,0 
EMUI)s0. 
FM(I)20,40 
5 P(IsJ)}=Q(134) 
el Le DE EE TT 
IF(KASES2) 35936935 
36 P(15»1)z1.40 
MEER 777 I 
P(353)=1-0 








-——_ _-_.___.__T—_———r ———————— —— S—— --- - a mm rm 
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PROGRAM OPTCON CONTINUED |. | | .— 


35 CONTINUE 
PRINT 19 
ΓΕ — 1595 FORMAT(U/72X5»56HNSTAGESGX»THAT (151) »3X» 7HAT(1»s2) 4X 56HP( 151)» — 
14Xs6HP( 192) »4X 5 6HP (2»51)»4X56HP(252)) 


e CEET E ο) 
DO . 22) KK=] s NS AGE E 
DEN=0.0 | 
DO6 I215»N 


6 EMCLI)SEM((IO)TDEL(J) *XP( JS I) 
DO8 ΙΞ1.Ν 
Eu" m DOT Jz DN € —— X ΠΜ 
7 FMCI)SFMOI) TEMCJ) *PHI(GCJS I) 
8 DEN=DEN+EM(I)®*DEL(I) 
DEN=-DEN=R 
DOLOISTSN 
AT(I)=FM(I)/DEN 
FM(I)=0%0 
C CALCULATE PSI(K»I»sJ) 
πο E 
DO13J215»N 
; 13 PSI(IsJ)}=PHI(Is4Y}+DEL(I)XAT(J}) 
e E PL... 
DO 15 IzloN 
DO 15 J=1sN 
DO 15 L*15N 
ο μυ τομ ο Όρο ο κο νο 
DO 17I=1sN 
PONI NEU  -""".——c— - 
DO 16 L=1sN 
16 HM( IsJ)=HM( Ts JI +GM(1»L)*PSI(L, J) 








— —— -- —-——— .- -——. 


C CASE 1 TERMINAL CONTROLS OMITSQOCT SJIJIIN EQUAT I ON “KORPUS 








C CASE2 MINIMUM FUEL OMIT Q(I»sJ) IN EQUATION FOR P(I,J) 
Enc OU pESRASE-2) 31»31333 
31 P(I3J)=HM(IsJ)}+RXAT(I)*AT(J) 
GO TO 17. 
c CASE THREE MINIMIZATION OF TIME AND FUEL 
33 PLISJ)SHMCISJ) QUI SJ) RKAT CI )&AT (J) 
hcec 7 9) S00 a μμ πω. ο 
DO18 Is1»N 
DO18 J=1>N 


ER Et Leet CR eg, ` n SS 
PRINT 20»KKsAT(1)»sAT(2)5»5P(1»1) PC 1»2) 5PC 251) 5P C252) 
20 FORMAT(I55(6F10.46)) | 
EE 22 CONTUR EE 
1111 CONTINUE 
END 


Ln Jen ————— 


ne eer ο 


SS sr va IT E EE = ——Á—————— (€ — e 





end e e e sn 


PROGRAM OPTCON CONTINUED ` 


SUBROUTINE PHIDELI(PHI sDELsNsDT) 
C VALID ONLY FOR A CONSTANT F MATRIX 
DIMENSION F(12512)5»PHI (12»12),TERM(12»12)5WORM(125»12) 
ls DEL(12) sDELM( 12512) »TELM(12,12)  DELP(125,12)»9D( 12) 
1 FORMAT ((8F10-0)) 
READISCCFUÜCIRSICIS,ICSI1SNOSIRS1,N) ` Le 
READ 1 (D(I1)s1=1,N) 
1003 PRINT 399sDTs(CFCIRs IC) sIC=loN) sIR=1 9N) 
ET SOT IJE /»((2F8,2))) 
PRINT 3991 (D(1)sI=lgN) 
3991 FORMAT(/ SHOtI)s/t2F' 8922) 
E A NFINALUSI 
TM=0 eO 
DO 400 IRz1,N 
DO 400 IC*15,5N 
TERMC(IRsIC)=0%0 
WORM(TRsIC)=0-0 
TERMCIRSIRYSISO 
TELM(TReIC)=TERM( IRs 1C)*DT 
DELPCIRSIC)STELMCIRSIC) 
e RO TT 
DEL (IR)20.40 
400 PHI(IRsIC)=TERMC(IRsIC) 
HE. e οσο 
DO 500 IR=l>N 
DO 500 ICzloN 
EM -— DO-500 "JNsT,N TTT ος 
DELM(IRsIC)=DELM(IRsIC)-TELMI(IRsUN)*F(UN»s IC) *DT/(TM+160) 
500 WORM(IRsIC)=TERM(IRs» ἽΝ ΕΟΝ» IC) *DT/TM*WORM(IR9IC) 
We POO] <i 5s. -  . -  -- °° SS eee 
DO 401 IC=1,N 
TERM(IR»IC)=WORMC(IRsSIC) 
MARN IREI C =DECMCIR = ——————— 
DELPLIRSIC) SDELPLIRSICIFTELMCIR SIC) 
PHICIRsIC)=sPHIC(IR,IC)+TERMC(IR, IC) 
DECIR; LE) =030 mm ` 
401 WORM(IRsIC)=0.0 > 
ABC=0.0 
ο Zen τε η ο σσ 
DO 2. Ξ5]1Ν 
AASTERMCIS9J) 
EN MENO) πα μη ME 
ΙΕ(ΑΒς-ΑΒ) 329292 
3 ABC=AB 
I NJE TT er "> io π΄ 
IRR 00 000005-ABC) 5,56 
_5 GO TO 4 
ο. G PRINT 11s((PHITIS ITS JEI NT; IFI, N — COCO 
LI FORMAT(/9X»9HPHI(I»J)=/(2F15611)) 
DO 600 Is1»N 


A MR s Hu a i E I E TEE MR RM 


—————— Á—— == —— e e 
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cun Tata 


- nio ti M a co — — G 


PROGRAM OPTCON CONTINUED 


DO 600 K=1sN 
DO 600 U=lsN 
5900 DEL(IJsDELOIOSPEUROGIESUNS XDEEPTJS, COLDEDINSS 
PRINT 12s Gye NI 
12 FORMAT(SXSTHOELETT)SZ(IF15.11)) 


END 
END 
η o ^- ——-—-— πο — ——  . — G ΤΑΞΗ a 
2 20 
1.0 
Ene l.0  — ^ m 2" — CC ^C" πα: 
1.0 =o eo 
360 
Em 6 
2 20 
1,0 Get? 
EE — — τ τω. Te 
1.0 -3.5 
3.0 
ELI uu ν᾽ 
2 20 
1.0 1,0 
ees "EE ae 
ito -3.9 
320 


tnt RI TT min E e ον 





r_1_1_rrTT_rrrEttttdtd _c _——e=-_—— e e Om oss cae 


Variable 


iS 


op 


DEFINITION OF TERMS FOR PROGRAM FILTER 


A 


Is Ia 


IN 


IN» Ix» 


Fei Χ- 


W(k) 


Definition 


nxn state transition matrix 
nxl control distribution vector 


nxl optimum filter gain matrix 


plant observability matrix 


measurement noise covariance 
matrix 


random excitation distribution 
covariance matrix 


error covariance matrix 


nxl vector of plant states 


nxl vector of plant non noisy 
observable states 


nxl vector of noisy plant ob- 
servable states 


nxl vector of the predicted 
value of X(k), based on the 
previous observation X(k-1) 


nxl vector of the predicted 
value of Z(k), based on the 
previous observation Z(k) 


nxl vector of the best estimate 
of X(k), based on the current 
observation of Z(k) 


random disturbance of environ- 
ment or excitation noise 


TABLE II-2 
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Remarks 


constant for a given 
time-invariant plant 


constant for a given 
time-invariant plant 


order of matrix depends 
upon number of plant 
observables 


value selected by 
operator 


value selected by 
operator 


initial value selected 
by operator 


assumed normal with 
mean zero 


Variable Definition Remarks 


V (k) random disturbance of assumed normal with 
measurement noise mean Zero 
Ow standard deviation of 


excitation noise 


σι standard deviation of 
measurement noise 


TABLE II-2(cont'd) 
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ROGRAM FI 
C DI ORDER OF SYSTEM IN I2 FORMAT 
imc DZ SAMPLING INTERVAL IN-8FIO.;O0 FORMAT 
È D3 F MATRIX BY ROWS IN 8F1040 FORMAT 
C 04 D MATRIX BY COLUMN IN 8F10,0 FORMAT 


τ D5 VARIANCE OF EXCITATION NOISE SIGWSQ IN 8F10,.,6 FORMAT = 
D6 R COVARIANCE MATRIX OF MEASUREMENT NOISE IN 8F104,6 FORMAT 
PHI SYSTEM TRANSITION MATRIX 
ENC  — DEL DISTRISUTTON MATRIX 
G OPTIMUM GAIN MATRIX 
H OBSERVABLE MATRIX t 
En P BEST ESTIMATE OF ERROR COVARIANCE MATRIX” = CI" ——— 7 
Q EXCITATION NOISE COVARIANCE MATRIX 
DIMENSION P(12912)»Q(12»12)»H(12912)»R(12912)»G(12912)»PHIT(12912) 
EE  I1,PHI(CI2*912)»DEL(OCI2)5DELDELT(12, 12) 5 DELSTI2»12) »DELST(I2, 12)» 
2X(12)sXHAT(12)sYHAT(12) »sPNEW(12,12) 
READ 339 N 
e 4121 Në 7 
READ 2,DT 
2 FORMAT (8F 10-20) 
ιο; c — E UR ------ 
1003 FORMAT(1H1) 
CALL PHIDEL (PHI »,DEL»NoDT) 
En DO IO] !ll*I1»4^ = = TTT 
READ 20»SIGWSQ 
20 FORMAT(F10.6) 
DO 1001 LK-1»4 — 
C INIALIZE MATRICES FOR EACH RUN 
DO 10 Izlo9N 
eee OJ SN ST EE 
PNEW(IsJ)=9.0 | 
P(I»$J)20,0 
ο κ) πο Ὁ —_— — 
Q(I,J)=0-0 
DELDELT(I»J)=0-0 
ορ Πο ο TT TT 
10 DELST(I»J)=0.0 
READ 200196R(191) 
2001 FORMAT (8F 10-6) | ε.α αμ. 
DO 30 I721,N 
30 DBleS(Isl)=DEL(T1I) 
CALU TRANS(DELS,N»N.DEUST) ~ oe τ 
CALL PRODI(DFLSIDELST NoNoN, DELDELT) 
DO 40 IzlonN | 
Ep DO 40"J*1,N ^ "Uc a a, m 
40 Q(I»J)35SIGWSQ*DELOELT(I»9J) 
SIGW=SORTF(SIGWSQ) 
E PRINT 1004,STGK TT 
1004 FORMAT(/ 10X»5HSIGWs /5E20.48) 


AND AA 


C 











LE Im CER EL .... 





-O— ————— ni ee > n “ ree  ———— - 
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a —— —— = ee — 


PROGRAM FILTER CONTINUED 


P(15,»1)2,02 
PLISZN SO 0 
P(2311*0.0 e ο πο... 
P(292)=1.-0 





PRINT 2002,;RCL1s1)5 (CQUISJ) sJ71, NO) » Iz 15 NOS CCPCISUJYS Je 1 ND » I T19ND 
0 — 2002 FORMAT(U77/720X38HRCISI)8 5F10.6/20X, 8HQTT5 7) 9 ,5FIO.67/20X5 BHPUCISJ] 
ΠΟ. 6) 
Η(1 91) 51. 
i HZ 
PRINT 700 
700 FORMAT(//5X53HG11» 7X» 3HG12» 7X, 3HG215» 7X, 3HG22» 7X » 3HP 11, 7X » 3HP 12 , 7X 
I SÜP2lIY.TX32HB2274y——————————— À» À c Wu ÁN 
DO 1000°KK= ier 
CALL GP(HsPHI »PsQsRo2919GsPNEW) 
PRINT—800CCS UC EE = TN = DN CP NEW Css = SN = ehre 
800 FORMAT (10F 10.5) 
pom EET. | 
— ——50—11—J- T5 lb ENN 
11 P(ISJ)SPNEW(CIS J) 
1000 CONTINUE 
OR e MM 
1005 FORMAT(1H1) 
1001 CONTINUE 
cip ee oc c Neo M E 


^ 


SUBROUT-INE—-GP-CHs PHI 9 P10 sROKN SK PPOs PNEW E NI 
DIMENSION H(12912) sPHI(12912),P(12912) sQ(12912)9R(12912) 96112912) 
IENEVWTIN 1298 HTQOI2s12).» TVLUI25 1202 B2 le) 
ολοι TRANS ΚΡΑΝΗ -- 
CALL PROD(PsHTsKNoKNsKPoTV1) 

CALL PRODIHSTVISKPSKNIKPSTVZ2) 

EE EE ΜΗ iii 
eE PKP n 0000000000001 TV E EERI 
IF(KER-2) 1015110, 101 | 

—1l10-PRINT—d1——— — -— : i ———————M——— 

111 FORMAT(5HKERz2) 

101 CALL PROD(HTsSTV2,KN5SKP KP, TV1) | 
CAEL—PROD (P^, TV, KN 9 KN 9 KP’ GJ) 
CALL PROD(H,PsKPsKNsKNsTV1) 

GAM DOBIIGSTVI  KNSKPSKNe TV2) | 
DO-r02-17=1, KN===" e - SS 
DO 102 Jz15KN 
ee πλω 
—CALI—ADDEPSTV25 KNIKNITVI )——— — 
ΓΕ ΠΙΤΤΥΙ,ΚΝ»ΚΝ»ΚΝ»ΤνΖΙ 
CALL TRANSIPHISKNSKNSTVI) 
— CAEL PRODIT V2; TV BRN KNGKNSEN Br — 
CALL ADD(PNEW+QsKNsKNsPNEW) 
END 





————————————— — 


M —X—9— CL À — ——— ———- 
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PROGRAM FILTER CONTINUED 


——————————————————— 








SUBROUTINE TRANS(AsNoMoC) 
DIMENSION A(125121;C(12;12) ° 
DO 153 I = 1)Ν 
DO 153 ER 19M 

EM 30153 C(J,I) = A(T8I) — 
END 


SUBROUTINE RECIPINSEPSASXIKER) 
DIMENSION A(12512)5X(125»12) 
lk EE ee — = ae 
DO l1 Js1»,N 
1] X(15J)20,0 
Eno — 3O—0DOZ2 KK ON 
2 X(KsK)=1-0 
10 DO 34 L=l aN 
En n. KPzÜ 
220.40 
DO 12 K=LsN : 
ΥΠ ΤΙΑΙΚΘΟΥΤΙΙΓΥΤΕν δ = = "eer" 
11 Z=ABSF(A(KosL)) 
KP = K 
12 CONTINUE 
ΠΕΤΕΡ 13920220 
13 DO 14 JzL,N 


— cc - © — 


=  ————Ùu o ro rn =. 


——— — — — —— — ——— 





ACLIJIZA(KPSJ) 
14 A(KP,J)2Z 
DORIS JEDN πι Ὁ- 
Z2X(L5J) 
X(LsJ)=X(KPsJ) 
qe sg —  —— 
Beer CABS UA(L»L))~EP )50550330 
30 IF (L~N) 31534934 
a 531 PS ELS] 
DO 36 K=LP1sN 
IF(A{KsL))32336332 
ee LEE eL)/ACLeaL) 
DO 33 Jetblah 
33 A(KsJ)=AlKsJ)-RATIO*®A(L »J) 
~=oilor_ 25>) — “TËS E 
35 XIKAOJ)ISXIKIJICRATIOKKELSJ) 
25 CONTINUE 
GON I NOT ~~ 
40 DO 43 1Ξ1.Ν 
[I=N+1-!1 
en ena 5 — | sN |... | <<... °° 
S=0.0 ì 








e -me e—— —- νο 


———————— 
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PROGRAM FILTER CONTINUED 


IF(II-N)41343343 

41 IIP1=II+1 

BE rech 

42 SZSTACIISK)EKXËKOJ) 

43 X(IIs4)=(X(IIs4)-S)ZA(IIsI1) 
Ee ee 
RETURN 

50 KER=2 
ENOT 0077 


SUBROUT INDa ge IDEE EHI DEN DT) 
C VALID ONLY FOR A CONSTANT F MATRIX 
DIMENSION F(12512)5PHI (12»12), TERM(12512) sWORM( 125912) 
E Ls DEE(C12T,DEDMU12* 12), TEEM(12,; 12) DECP C12, TONO" ΤΠ ΝΙΝ 
READ1s((F(IRsIC)sIC=1s0N)s IR=1gN) 
] FORMAT ((8F106%0)) 
READ (Dre SN TIT 
1003 PRINT 399,0T$CUFCIRS IC) »ICS1sN) »IRSISN) | 
EXSCOUOERORMATDUAZZSmo»Dzs.1Eb.3///.57HF(I,,J)s7,00t8F8v2)99 
PORE O IeLSNT-———— —ÉÉÉÉE S 
SEE OAM DI  BBFSa2)) 
NFINAL=1 
Ee ___ 
DO 400 IR=1,N 
DO 400 IC=1,N 
Εκ μι κ 1.) =O60 - 
VORMCIRSIC) 0.0 
TERMCIRSIR)Z1.0 
ΙΤ ΠΠΡΓΚΗΙΓΟΤΧΡΤ--τ; 
ο ο πο με κ, [ο 
DELM(IRsIC)=0%0 
Ee U πο στ -u 
400 PHICIRSIC)ZTERMCIRSIC) 
4 TM=1e0+TM 
BOS 00 Ek Nm Ra Ene 9 
ΡΟ. 500 [IC=15sN 
DO 500 JN215N 
DEEMCEFRSTC)ZDELMUIR, IC)-TEEMCUIRS,JN)*FCAJNS IC)*DT/Z(UCTM* TO) 
500 WORM(IRsIC)=TERM(IRsUN)I*F(UNsIC)*DT/TM+WORM(IR9sIC) 
DO 401 IR=1,N 
ae Cer M as ii 
TERM(IR»sIC)=WORM(IRs IC) 
TELM(IRsIC)=DELM(IR,IC) 

OE DES PRIT E NR FF EB μμ | 
PHICIRSICOZPHI(CIR» EL IC) i 
DELM(IRsIC)=0+0 | | 

qUI- WORNWEEXGIRMEOCERU.O" X cuc | EN 








ABC=0.0 
DO 21z19N 


——————— —————— — — À ———_s — 
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PROGRAM FILTER CONTINUED 


DO 2J21»N 
AAZTERMCISJ) 
AB=ABSF (AA) TT, — — 
IF (ABC-AB)23,352 
_ 3 ABC=AB 
νην... 7 CONTINUE 0 — 9mm 
IF(0,000000005-ABC)55,556 
5 GO TO 4 | 
SeePRI ND DOZe((PRICIRe 1 GiaaloG@ iam. !R=lsNy TT E 
502 FORMAT(////9X»8HPHI(IsJ)///(6E2048)) 
DO 600 Is1,N 
i DO 600 K-*15N 
DO 600 «1.Ν 
600 DEL(I)sSDEL(I)TPHI(IS J) *DELP(J, K) *D(K) 
eee PRINT 503 (DEL(I),ISITSN)" "^  — 
503 FORMAT(////9Xs6HDEL(1)//(6E2068)//) 
END 


SUBROUTINE ADD (ASBSNSMIC) 
in — —DIMENSTON A€(12»12)»B(125»12)5C(12»12) 
DO 152 Iz21»N 
DO 152 JU=lsM 
e 2 CI) = AlCleJ) + Ες νο == TTT 
END 


———— — -—- - —- -———— M — —— —  —— — - 


SUBROUTINE PROD (AsBoNoMoloC) 


DIMENSION A(12912)58(12912)9C( 12512) 
bose, jj et rn 


DO 151 Jz15L 

Cie? ΞΟ. 

DO 151K - 1s5M 
151 C(I»J) = C(ISJ) * ACISKO*B(OK9J) 





————— en οι EE Pl D G A e E παπα παντα, 





END 
gi a m 
2 
E — — — — 
1.0 -26.5 
360 
M G OO e Tr —— 
«0 
202 
κ] E — €— ᾽ς -- 
1,0 


——— - ———— SIS TI Su ---- -ᾱ- --- -- 
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c PROGRAM HYBRID II 


CC —TH1S-PROGRAM..USES VARIABLE FILTER GAINS e 
C OTHERWISE IT IS FSSENTIALLY PROGRAM HYBRID 1 
i 






———— — — IHILS PROGRAM CLOSES THE-.LOOP--ON-THE--OPT-LMUM-F-IL T-ERCONTROL-L.ER—— — 
C PROBLEM. IT ASSUMES THAT STABLE CONTROLLER AND FILTER GAIN MATRICE 
C HAVE BEEN COMPUTED ON THE BASIS OF DESIRED RESPONSE AND THE 

CC STATISTICAL PROPERTIES OF. THE ANTICIPATED RANDOM DISTURBANCE-ANO— 





Ge MEASUREMENT NOISE. THE SYSTEM MAY BE DRIVEN EITHER BY INITIAL 
δ CONDITIONS OR BY STEP OR RAMP INPUTS,OR BY ANY COMBINATION OF 
E E A o_o μμ. —- eme 
e THE PROGRAM SOLVES THE FOLLOWING | EQUATIONS 
C Y(K)sH*X(K) 
mecc UN Ss La AMET Ketter omg E τ 
a XS(K}=(I*GH)IPHI*XS({K>1)+(1I= -6H) *DEL*ATUXS(K- 1)-DINP(K-71)) *G*Z(K) 
C X(K+1)=PHI*X(K}14DEL*W(K)+DEL#AT*X(XS(K)=DINP(K)) 
— C... WHERE AM 15 MEASUREMENT_NOISEsW IS THE RANDOM DISTURBANCE SAND — 
C DINP IS THE DETERMINISTIC INPUT | 


DIMENSION X(12912) 9XS$(12912)sSIGV(12)9Y(12912) 9Z(12512) » PHI (12512) 
a DEL (12512) 5H(612512) sAT(12 512) 96(12912) 3 TEMP1(125 12)  TEMP2(125.12)5—4 
2TEMP3( 12912) sTEMP4(12512) sTELP( 12912) sTELP1(12512),TELP2(12512) > 
3V(12912) sDINP( 12912) 91 T(12)9TMF (900) 9X1(900) »XS1( 900) 9X2(900), 
4X§2(900)5U(12512),B(12912) »P (12912) δι 129112) 81119000) 1900)  — 
READ 1s(IT(I)»,1=1,6) 
READ εί Γον E) 
δω doi FAV NE TS n ^ — --- 
DO 1111 Ks1,5 
6 INPUT SOME CONSTANTS» AS NOTED 
í— Tele. -—" È "a a 
e KN IS SYSTEM ORDER. KP 159 THE NUMBER OF OBSERVABLES. 
KN= 2 
PIE= 3. 1415926 
AR=0 5 
rr THËTARSPTËZE RO — | E e 
APSOri 
THETAPSPIE/255 
etn eee’: | te Pent Ge MAGN TILDE OF. STEPSRAMB. LÉ£ MAGNITUDE DE RARR:a RE Laser πμ... 














C AND PITCH SIGNALS ARE OUTPUT OF ZERO ORDER HOLD, CONSTANT OVER SAMI 
STEP=4,0 
cette AMET nn | 
C DT IS INVOLVED IN UPDATING ROLL OR PITCH INPUTS - 
DT=0%0 
— OR, IS INVOLVED IN UPDATING RAMP ÍNPUT — 
DR=0.0 
C DA IS USED TO REDUCE STEADY STATE ERROR TO RAMP INPUT TO ZERO 
cu —— POR. LHILS-LYPE-QNE. BYE TEM ii, gine 
DA=0.0 


E SETTING DINP= "m 0 AT THIS POINT INSURES THAT NO DETERMINISTIC 
€ INPUTS EXIST PRIOR TO TIME = ZERO 

(NANA I IN 

DINP(231)=0%0 
AG INITIALIZE MATRICES 

EE, LS kal __________o oo __ 
DO 66 J=15KN 
Y(1IsJ)=0-0 
Z-( T-sJ-)-=0 4. O- 
UtI»J)2040 
B(I»J)20,0 
Τ᾿. ΓΞ ο ο ως 
G(I»J)=0.0 
H(I,5J)20,0 

o c o  --. 
66 PHI(I»J)2040 

DO 65 I=1,KN 

MENU EE STOVU 0e 0  ——————. 
REAO 2»SIGV(1)5SIGW 

2 FORMAT(2F 10.6) 

e logi) = Se NITTO), — — — sa en 

SIGWSQzSIGW*SIGW 





C CARDS ABOVE PROVIDE FOR SETTING DESIRED VARIANCES IN W AND Vo 
EC  —S$IGW—L$—THE-STANDARD-—DEV- A.-T-ON—OF-—W .—S I M-CARULY FOR—SIGV 

C 

C INPUT MATRICES 


dee O a mumu ο”... 
PHI(15,2)204,277086 
PHI(2»1)2040 | 
Mn etal er. 7 JamÜe03019J.— __ ——— M i ---- EE 
DEL(131)=0+619640 
DEL(2»1)204831259 
I 1 . 1L —————M———— τ. 
H(1»2)20.0 
P(1»1)520,02 
οι L0.0— ——— Bed Toa E ee 
P(2,1)=0600 
P(252)21.40 
AT(191) 2-14-2030 : Si 
AT(192)=-0,3426 
C 
CE THE. FOLLOWING CARD INITIALIZES THE RANDOM NUMBER GENERATOR. .. .. .— 
NUNIF=1220703125 
PRINT 3 
ἔν... 3A Li L/ 
PRINT 800 
800 FORMAT(//4Xs4HX(1)»4X»5HXS(1)»5 3X» 4HX(2) »54X»*5HX5(2)5 3X» HV» 7X 3 1 HW, 
1-7X.»-3HG61-1-3 5X» 3HG12» 5X »4HAT11» 4X5 4HAT12»4X» HU, 7X 93 3HB119 5X, 3HB21»2X,- 
26HOINPi115,2X,6HOINP21 //) 


88. 


Oem o e ZATTERA Dm ORNAGO eer 


E SH 


meg OV ποιον μα. _ _ a l i 


g 


PROGRAM HYBRID II CONTINUED 


PLANT INITIAL CONDITIONS 


κ» »1)20 e). —— ee ee eee 


X(151)=0.20 

FILTER INILLAPENSONDITIONS 
CALL PRODIH»XsKPoKNo KP »¥) 
DO 12 I21,KP 

CALL RNDEVINUNIF sDEV) 


CALL ADD(YsVsKP919Z) 
Ἀθι η σε eee 


cn E 280, pa ο 


pee ODO 500 KK aie ee 


LO 


B- der ΕΙ; a o o 


1 


LL=70 





CALL FILTER(KN4sKPsPHI DEL sSIGWSQsRsHsPsG) 
CALL PROD(HsXsKPsKNslsY) 

DO 10 1z19KP 

CALL RNDEVINUNIF DEV) 

V(Isl)=SIGV(I)*DEV 


CALL ADD (etek Pol 7) o RI OE e TSE 


CALL PROD( PHI sXSsKNoKNols TEMPI) 
CALL PROD(GsHsKNsKPsKNsTEMP2) 





DO 11 Jz1,KN 

TEMPZi PJ) 

CALL PRODITEMP2,$TEMP1sKNsKNsl1sTEMP3) 
CALL ADD(TEMP1,TEMP3,KNs1sTEMP3) 

CALL ADD(XSsDINPsKNs1sTELP) 

CALL PROD(ATsTELPs1lsKNs1lsTEMP1) — 
ΤΈΜΝΕΙ ΠΤΕΜΡΤ1(1»1) ΠΑ 

CALL Ὁ Go PELSTËEMP1SKNI1: 15 TEREN 


ΓΑ. ΡΚΩΡΙΤΕ5ΡΘΤΕΙΡωΚΝ»ΚΝο1«ΤειρΡ1], --᾽---- 


CALL ADD(TELPsTELP1sKNs1sTELP1) 

CALL PRODI(GsZsKNsKPs1sTELP2) 

CALL ADDITEMP3STELPLËSKNS19XS ) - ἃ 
CALL ADD(XSsTELP2sKNs19XS) 


ΓΕ e SO KR °° dl. 


IF(KK-4) 17517518 


18 CONTINUE 


—— — rr——— 


-æ 


" RAMP=1.-0 _— n 3 — ee 
DRSDRY10 
DASRAMPK(36 573.0) 
NE I 
DINP(1» 1)=- (STEP+RAMP*DR+AR*SINF(THETAR*DT)) 


DINP(291)=-RAMP 


ae  . 





XI(KK)=X(191) 
ASS ll) 


— — - — —— ss ου — oo — 
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PROGRAM HYBRID II CONTINUED ua. 


X2(KK)2X(25»1) 
XS2(KK)=XS( 291) 
B11(KK)=Bl(151). .. 
UU(KK)=U( 151) È 
EDENICSOLSXCLsSI)9S(1»1)5X(2»199XS9S(291)9VUI21)9»W»GUls)»G(.29M) s 
SS SSSLLAT.( 10.1.) s AT(192) pUl(.10.1.).95( 191.) 0B(201)0DINP(141)-0DINP (2.01) 
801 FORMAT(15F8,4) 
CALL RNDEV((NUNIF ,DEV) 
MeSIGWWDbEV ______—_——_—mm_ 
CALL PROD(PHI sXsKNoKNolsTEMP] ) 
DO 803 Iz1»KN 
we 803 TEMP2(151)=WeDEL (1 51) 
CALL ADD(XSsDINPsKNs1sTELP) 
DO 13 ΙΞ1 ΚΝ 
T DO 13 JE ËN .. 
Pe Bi led) =TELPC I oJ) 
CALL PROD(ATsTELPslsKNolsTELP} ) 
DO 14 I21».KN .. 
DO 14 Jz1,KN 
14 U(IsJ)=TELP1(134) 
= — [ELPI(ls1)sTELP1(CIIA1)9DA ἃ νο”... 
CALL PROD(DEL»sTELP1lsKNolols TELP2) 
CALL ADD(TEMP1sTEMP2 sKNol 5X) 
CALL ADD(XsTELP2sKNolsX) |. | NT 
900 CONTINUE 
MC=1 
a Lc PN eo --. 
CALL DRAN ILL »TME» X1 sMCs0sLAsIT9100091009390929290910091sLAST) 
iC=2 
LA=4HXS1].. | |. |. m == es ATE 
CALL DRAW(LLs» ΤΗΣ χο Με ο aO gë toa ume 0»2»2»09»109»1» LAST) 
MC=2 
LAZHHX2 ..— I | 
CALL DRAW(LL» mea e Më ere is 06093» "UTR 10919LAST) 
MC=2 
ee LA =SGHXS2 m è i IA 
CALL DRAW(LL »TME» XS2 MC 9 Gre Ee, 090% on 3% PIERO.) 
MC=3 
LA=4H31]] | | |). slot? DL το - 
CALL DRAW( LL» TME, Blw ῃ EE, 090% On e OO 1 Ga ASH 
MC=1 
' LAs4HX1-T _ JM xL 
CALL DRAW(LL» ΤΜΕ»Χ1.΄ EU on nu lOS oh MOS oo rosi LASTI 
MC=2 
LAs4HU-T ο απ ο ιο. 
CALL DRAN ILL TMESUU Ποπ... στη Ρα. 
MC=2 
LA=4HX2=T_. | cc cx. x . c 
CALL DRAW(LL»TME> X2 ` ANE EE an ο ο 2i12450 » lOsl»LAST) 
MC=3 


= PP——€—€—É«+—»____ - -—— 


fe. -e emea m me e a e o e a e mem ee re 








e a — - ege 








—— ut  -— 








PROGRAM HYBRID II CONTINUED zu 


— — — — ———- — — —M—M 


Te 4HX1X2 

CALL DRAW(LL 9X1 »X2 9MCSOSLASITs1eO 9 1e0O 9390 252 90910919LAST) 
ο CONSINUE e A 

END 


E — —SUuenOUT-ING-—PROD-—CANMB ήν. ο —— — 
DIMENSTON AUCY299995 558012512) 5,5C(12512) ^ 
DO 150 90-50! 
EMEN D05. 151. 9-15, 
C(IsJ)=06 i 
DO o NR UM 
po oe E — np - -- ϱ 
END 





e __SUBROUTINEC ADD (As B4N 5M C)_—___ o or to i 
DIMENSION e E eer ( 12912) 
DO 152 ΙΞ]1.Ν 
Ες ο πι E, 
L527. CCI 9pJy) AU EOM ES UD) 
END l 


SUBROUTINE FILTERI(NsKPsPHIsDEL ,SIGWSQsRsHsPsG) 
— C PHI SYSTEM TRANSITION=MATRIX — ~ 
C DEL DISTRIBUTION MATRIX 
C G OPTIMUM GAIN MATRIX 
ννν. ιν... OR SER Ri. o ë- ë ü ë 
C P BEST ESTIMATE OF ERROR COVARIANCE MATRIX | 
C Q EXCITATION NOISE COVARIANCE MATRIX I 
x — DIMENSION BR4(1-2»123)»Q:-12.9».112.).9H (12:91-29:12 911 2-)9 6123-12)» 2 ELI TTD πμ 
ΙΡ »129.,DEI«012)9DELDELT(12512) 9 DELES4 12912)» DEI SITIIZ, 1:208 
2PNEW(12512) 
E 1 7, 
30 DELS(IsS1)-2DEL(I) 
CALL TRANS(IDELSS,N, NS DELST) 
CALL PROD(DELS»sDELST»sNsNoNsDELDELT.)—.-- 
DO 40 Is1,N 
DOM ONU IN 
40 Q(IsJ)=SIGWSQ*DELDELT(.134) 
νι. ρήτρα, Ν.ΚΡ,α ΡΝΕΙΠΙ 
DOJ I qe 
E | mm 
11 PLISJ)SPNEMCISJ) 
END 





SUBROUTINE GPIUHISPHISPAQIRIKNIKPSGIPNEV) 

e  — _ __DIMENSFPON Hl124Ì 2) PHI (12512 P.(_1 291-220) QU 2 5-1-2) - RC 2s G1] 22M 
LPNEW( 12912) sHT(12912)9TV1(12912) sTV2(12512) 
CALL TRANS(HsKPsKNsHT) 


——————————— t 





PROGRAM HYBRID II CONTINUED 


CALL PRODIPSHTIKNSKNIKPSTVI) 
CALL PROD( Hs TV1lsKPsKNoKPoTV2) 

E ee oe) mamoa#î”" 
CALL RETCIP(KPs 0000000000001 »TVlsTV2sKER) 
IF(KER-2) 10151105101 

EA E ellen dl _rrrrr lg —— 
111 FORMAT (5HKER=2) 
101 CALL PROD(HT sTV2sKNy KP sKPoTV1) 
CALL-PRODIPSTVIOKNIKNIKPOIGJ— 
CALL PRODI(HsPsKPsKNsKNsTV1) 
CALL PROD(G,TVI1,,KN5,KPS,KN»TV2) 
EE 5S3 DO 1021 saga N 
DO 102 J4=1sKN 
102 TV2(134)==TV2(134) 
e .. ...........- 
CALL PRODIPHISTV1SKNISKNIKNITV2) 
CALL TRANS(IPHISKNISKNSTVI) 
CALL PRODLTVASTVIS KNË:SKNISKNIPNEVJ).. — 
CALL ADD(PNEWsQ,KN,KNsPNEW) 
END 
SUBROUTINE TRANS (AsNoMoC) 
DIMENSION A(12512)5»C(12»12) 

















νιν... Mee a N .............. 


DO 153 Y= lo 
153 CORSI)  AXIs.J) 
ν.. cp -T € —e«—  — .. EM. m 


SUBROUTINE RECIPINSEPSASXIKER) 


Enc — 31M. DEMENSION.AtC12512)52X(129122.  — —_ —_ _ 


DO 1 I=lsN 
DO 1 J4=1lsN | 
os __ MM — -- 
i DO 2 ΚΞ].Ν 
2 X(KsK)=1-0 


MEN | Ore ee ee — À——— ——————— À——————— 


KP=0 
Z=0%0 
EC PO KEL SN Ἵν ο lo Jo ....- Ee 
Ze SA lim) ) ) 1191242 
11 ZSABSFCACKSL)) 


EEN d ee E eee ee. 


12 CONTINUE 
PEUE-KP)15,20w20 


OS ge me 5 5 2. 


Z2A(L5J) 
A(LsJ)=A(KPsJ) 


E 14. A(KP.J)52. — . n a -o — 


DO 15 Jz1,N 
Z2X(L»J) 


e x VENI)  __ _ _ __o__O_____u ur _r 


= Sail amm 


PROGRAM HYBRID II CONTINUED 


15 X(KP,J)2Z 
20 IF(ABSF(ACLSL)) -EP)505 50, 30 
30—LF-CU-N )-31-5 34-934. — —— — 
31 LP1=L+1 
ΡΟ 36 ΚΞ“ 
dë (Ai E ai 
B2 RATIOZA (K ZEN LO 
DO 33 EINE 
33 — A-(-K-9-J-)-— A-(CK- J-)-7— R-A-T-I-O3-A-CL.-9 J-) ————— 
DO 35 Jz1,N 
35 X(K»J) SXO(OK Ss JJ SRATIOXX(UL 9$ J) 
36—CONTINUE 
34 CONTINUE 
20-—DO 4299 00 
NS e eee EEN 
δο CA ON 
S=0 250 
LEGLIZ-NO)41.,43,43-—— EE dle 
41 IIP1=II+1 
DONA 2 KOPJEN i 
SË — 
TOS SOSIO 
KER=1 








NE RES SN Sa CADE TERMINISTIC ΙΝ’ ΜΗ 


I ARA DE LE Gal Nd eh, (bile e 
e1414 0.1414 
elëlë 0.228 
de ne ο καν 0, 350 
el4l4 Oe D2 
01414 Oe 722 


eee — EE =o — 


—Q— E 





APPENDIX III 
PROGRAM OPTIMUM CONTROLLER 
III - 1 General Description 

Program Optimum Controller is a CDC 160 digital computer 
program written for the control of the simulated hybrid system 
discussed in Section 11 and shown in block diagram form in 
Figure 11-1. 

This program computes a control for each sampled interval 
given by 

u&) - AT [xe - pro9 | (111-1) 
where A! is the feedback gain matrix determined by the Fortran 
program 'Optcon' and is manually entered. 

X(K) is column vector describing the states of the plant 
that are sampled by the ADC once each sample interval, 

DI(K) is column vector describing the deterministic inputs 
that are generated on the analog computer and are sampled by the 
ADC once each sample interval. 

The program considers all matrices as square matrices and 
will handle up to fourth order systems as it stands, and higher 
order systems with an increase in available storage cells for 
the variables and temporaries. 

The program uses subroutines delay, mulop, matadd, neg- 
mat and matmul which are 12 bit variable decimal point routines 
that are described in Appendix VI. 

III - 2 Overflow Provisions 

Both product and summation overflow error halts are pro- 
vided in the matrix subroutines. There is no provision for over- 
riding an overflow, however, the scale factor can be changed in 
cell 27 to allow larger numbers to be used. See signal condition- 
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ing, Section 9, for further information on number scaling, 
III-3 Internal Timing Control | 

The setting of the internal timing control is explained in 
subroutine delay in Appendix VI. For the second order plant 
and a time delay of one second, a coarse setting of (0012), 
and a fine setting of (0000), worked reasonably well. However, 
better accuracy can be obtained By using external timing control. 
III-4 External Timing Control 

Figure III-1 shows the operation of the external clock with 
the digital computer and the ADC unit. The input disable jack 
of the ADC unit is connacted by means of an AND-gate to the 
input ready line of the CDC 160 computer input cable and pro- | 
vides a means of delaying the CDC 160 computer by an external 
timing device. To input a digital word from the ADC, the CDC 
160 computer sends an input request to the ADC at which time 
the ADC converts the analog voltage to a digital number on the 
selected A/D ομβήῃει. If, however, the input που. is heldat 
. ground level by means of an external m the input ready 
signal to the computer is delayed until the external device 
drops the input mee from ground to -3 volts. Upon com- 
pletion of the conversion process, the ADC sends an input 
ready signal to the computer after which the computer will in- 
put the digital word. 

The — clock puts out a pulse at the start of each 
sample. The length of this pulse can be varied, thereby chang- 
ing the length of the enabling window for the ADC. This is 
done by changing the value of a capacitor which controls the 
switching time of a monostable multivibrator. As a guide a 
0.1 micro-farad capacitor produces about a 280 micro-second 
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pulse width. During each sample cycle, as shown in Figure 
III-2, the ADC must be enabled long enough to allow all 
channels to be sampled, and inhibited for the remainder of 
the cycle. Then the computer performs calculations to obtain 
the desired control, outputs this control D/A, and is hung up 
waiting for the next clock pulse. 

The enabling window must also be short enough so that 
the next set of samples is not taken until the sample interval 
has elapsed. The internal time can be adjusted to protect 
against this possibility by using subroutine delay. For this 
second order plant, . the coarse control was set to (0001), 
and the fine control to (0000)g. 

Π1-5 Use of Program Optimum Controller 

A. Set up the plant to be simulated on the PACE TR20 
analog computer similar to Figure 11-1 so that each state is 
sampled by the ADC. Connect the three remote control con- 
nections of the analog computer to those of the ADC so that 
the master clear switch on the CDC 160 digital computer has 
control of the hybrid simulation. 

B. Turn on the ADC/DAC, and the CDC 168 arithmetic 
unit. 

C. Load the program at (0000) and the subroutines at 
addresses specified by the directory cells. 

D. Enter the constants that change with the plant and 
the order of the plant as designated by an asterisk in the 
following listing of the program. 

E. Set the delay loop or the external clock for the sample 
interval desired, and run the program from (0000), . The 
following listing is for the second order gun train plant of 
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EXTERNAL CLOCK CONTROL 





INPUT DISAELE JACK 








Y 
-2y ADC_INIT LOGIC LEVELS 
| ! AND GROUND È ο: 
| CDC 160. INPUT REQUEST | INVERTER L CATE 
i | CINE mer -2V + 
| COMPUTER [DATA READY | | 
: LINE 





A/D CONVERSION 


FIGURE IlL-L 


TIME SEQUENCE OF ONE SAMPLE 


A/D COMPUTATION D/A SUBROUTINE 
ENABLING OF OUTPUT OF 
WINDOW CONTROL CONTROL DELAY 

| V | | 






VOLTS 






% WIDTH OF ENABLING WINDOW APPRO 
280 usec for each Oel uf CAPACITOR 


I— ONZ SAMPLE CYCLE — 


FIGURE LI 
SE 


-3 





Figure 4-2 vvith a sample interval of one second provided by 
an external clock. Three 0.1 micro-farad capacitors in 
parallel were required to make the width of the enabling win- 


dow for A/D sampling long enough. 


IS 


PROGRAM OPTIMUM CONTROLLER 





CELL MACHINE SYMBOLIC REMARKS 
CODE 

0000 ZG? jfi l Run from zero 

0001 1000 1000 location of program 
constants to be loaded as 
noted 

0017 0020 counter for null matrix 

0024 0004 *product m*n 

0026 0000 number of stages, set equal 
to zero for continuous 
sampling 

0027 0400 scale factor for mulop, 
scale factor of three 

0030 0002 *m matrix dimensions 

0031 0002 *n matrix dimensions 

0032 0002 *p matrix dimensions | 

0037 0000 master counter for number 
of samples 
addresses of subroutines 

0052 4300 delay 

0053 4400 mulop 

0055 4600 Ἴ matadd 

0056 4700 Ww negmat 

0057 5000 natmul 
addresses of variables 
and temporaries 

0100 to 0117 null matrix 

0120 0000 X(1,1) plant position 

0121 0000 

0122 0000 X(2,1) plant velocity 

0123 0000 

0140 0000 DI(1,1) position input 

0141 0000 

0142 0000 DI(2,1) velocity input 

0143 0000 


0160 to 0177 
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-DI (K) 





PROGRAM OPTIMUM CONTROLLER CONTINUED 
MACHINE 


CELL 


0200 to 


0220 


0221 


0222 
0223 


0240 
0241 
0242 
0243 


CODE 


0217 


7313 


7650 


0000 
0000 


0000 
0000 
0000 
0000 


loopl 


SYMBOLIC REMARKS 
templ 
*AT (1 ,1) position feedback 
gain 
*AI (1,2) velocity feedback 
gain 
U(K) scalar control 
ldc null matrix equal to Zero 
0100 
stf a 
lede. l 7 -20 
std 77 
ldn 0 
stm 
XXXX 
aob a increment address 
aod 77 increment counter 
nab loopl 
pta initialize X(K) = 0 
jpi matmul 
oloo 
0120 
0120 
pta initialize DI (K) = 0 
jpi matmul 
0100 
0140 
0140 
pta initialize U(K) = 0 
jpi matmul 
0100 
0240 
0240 


95 


PROGRAM OPTIMUM CONTROLLER CONTINUED 


CELL 
1032 
1033 


1034 
1035 
1036 
1037 
1040 
1041 
1042 
1043 
1044 
1045 
1046 
1047 
1050 
DUST 


1052 
10.5 
1094 


1055: 


1056 


1057 
1060 
1061 
1062 
1063 
1064 
1065 
1066 
1067 
1070 
1071 
1072 
1073 


“ra 
br 


MACHINE 


CODE 


2426 


4037 


2430 


4077 
2200 
1402 
4205 
2200 
0120 
4205 
7500 
XXXX 
7600 
4100 
XXXX 
5704 


2030 
5303 
5477 
6511 


2430 


4077 
2200 
1407 
4205 
2200 
0140 
4205 
7500 
XXXX 
7600 
4100 
XXXX 
9704 


loop 


loop2 
b 


C 


loop3 
d 


e 


SYMBOLIC 
lcd 26 
std 37 


lcd 30 
std 77 
ldc 

1402 
stf b 
ldc 

0120 
stf c 
exf 

XXXX 
ina 
stm 

XXXX 
aob b 


ldd 30 
rab c 

aod 77 
nzb loop2 


led 30 


std 77 
ldc 
1407 
sti d 
ldc 
0140 
stf e 
exf 
XXXX 
ina 
stm 
XXXX 
aob d 
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REMARKS 
set cell 26 = 0 for continuous 
sampling 


master counter 


input X(K) plant states 


analog to digital channel 


address of X(K) 


select ADC 


increment channel, address 
& counter 


input DI(K) deterministic 
inputs 


analog to digital channel 


address: of. DI (K) 


select ADC 


increment channel, address 
& counter 








PROGRAM OPTIMUM CONTROLLER CONTINUED 


CELL 


1074 
1075 
1076 
1077 


1100 
1101 
1102 
1103 
1104 
1105 
1106 
1107 
1110 


1111 
Tr 
μις 
1114 
1119 


1116 
1117 
1120 
ΜΙ 21 


el 22 
1123 


A mm ci 


MACHINE 
CODE 


2030 
9303 
5477 
6511 


0101 
7056 
0140 
0160 
0101 
7055 
0160 
0120 
0200 


0101 
7057 
0220 
0200 
0240 


7500 
2411 
79) 
0241 


0101 
7052 
0012 
0000 
5437 


6573 
7700 


0240 1 


SYMBOLIC 


ldd 30 
rab c 

aod 77 
nzb loop3 


pta 


jpi negmat 
DI(K) 
-DI(K) 


pta 


joi matadd 
-DI(K) 
X(K) 
templ 


pia 


jpi matmul 
AT 
templ 
UK) 


exf 


2411 
out f 
0241 


pta 


jpi delay 
coarse 


fine 


aod 
nzb 
hlt 


` i Sn o o_o ri; — eem 


REMARKS 


DI (K) = -DI(K) 


templ = X(K) - DI(K) 
U (K) = al EG - Di (9 |seatar 


output U(K) digital to analog 
channel one 


last word output address 
plus one 


internal time delay 

for external clock control 
set 

coarse control equal to 
(0001)g 


increment master counter 


first word output address 


APPENDIX IV 
PROGRAM HYBRID I 
IV-1 General Description 

Program Hybrid I is a CDC 160 digital computer program 
written for the control of the simulated hybrid system discussed 
in Section 16, shown in block diagram form in Figure 16-1, 
and in flow chart form in Figure IV-1. 

The program uses a filter to estimate all the states of a 
plant and a controller to synthesize the desired control. The 
output of the program is a control for the plant simulated on 
the PACE TR20 analog computer. Stable gain matrices are used 
for the filter and the controller in Hybrid I. The program con- 
siders all matrices as square matrices and will handle up to 
fourth order systems as it stands, and higher order systems 
with an increase in available storage cells for the variables 
and temporaries. The program can be used for plants with 
more than one observable state. The following listing of the 
program is for the second order train plant of Figure 4-2 with 
a sample interval of one second and with à pseudo noise-to- 
signal power ratio of one taking the noise R(1,1) = 0.02. An 
enabling window width of 630 microseconds is sufficient to 
allow all the A/D sampling to be completed on each pass. 

The program uses subroutines, xstar, listop, delay, mulop, 
adop, matadd, negmat and matmul which are 12 bit variable 
decimal point routines that are described in Appendix VI. 

For overflow provisions, internal timing control and ex- 
ternal timing control see Sections 2, 3, and 4 of Appendix III 
for program Optimum Controller. 
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Enter constants, program and subroutines 
STEP DA =0 enter A with 0000 
RAMP DA KO enter A with 4000 (for type one plant 


Ν/ 


EN 


Initielization of metrices 
Update G & P = PNEW for Hybrid II only 





initialization of Noise 
YES enter A with 0000. 
NO nitor À with OO 


V 


Turi 





τ 


Input ADC Y(K) ana DI(K) 











Calculate V(K) = sigvimoise. 
| Z(K) = VOR) + ¥(K) 
W(K) = sigwinoise 


stage onim v 
me | DIS AS(K) 
———— 
| 
| 
| 
| 
| 


U Galcalate U(K) =A? [ks(x)-DI(K] — Seet) 
Output UÇK) 


| 







halt 







| Update matrices: DI(K-1) = DI(K) 
; G and P = PNEW for Hybrid II only 
p—————— 











internal or External 


FLOWCHART FOR HYBRID PROGRAMS 
FIGURE IV-l 


- h 


IV-2 Use"obbuscmubybrid- 
A. Set up the plant to be simulated on the PACE TR20 analog 


computer similar to Figure 16-1. Connect the three remote con- 
trol connections of the analog computer to those of the ADC so 
that the master clear switch on the CDC 160 digital computer 
has control of the hybrid simulation. 

B. Turn on the ADC/DAC unit, and the CDC 168 arithmetic 
unit. 

C. Load the program at (0000) and the subroutines at 
the addresses specified by the directory cells. 

.D. iEnter the constants that change with the plant and 
the order of the plant as designated by an asterisk in the 
following listing of the program. 

E. Set the ES loop or the external clock for the sample 
interval desired. | 

| F. Since the gun train plant is a type one system it has 

a steady state error when responding to a ramp input. A DA 
factor was introduced to nullify this steady state error. Hence 
for ramp inputs enter A with (4000), run the program from 
(0000). and halt in cell 1047. Enter A with (0001), to not 
reinitialize the noise table otherwise it is restored to the same 


order, and run the program again. 
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PROGRAM HYBRIDI 


CELL MACHINE SYMBOLIC REMARKS 
CODE 

0000 7101 jfi 1 run from zero 

0001 1000 1000 location of program 
constants to be loaded as noted 

0017 0020 counter for null matrix 

0021 0072 *standard deviation of environ- 
ment noise 'sigw' 

0022 0044 *standard deviation of measure- 
ment noise 'sigv' 

0023 0044 *number of observables 

0024 0004 *produce m*n 

0025 0000 *DA steady state error correction 

0026 0000 il number of stages, set equal 
to zero for continuous sampling 

0027 0400 scale factor for mulop, scale 
factor of three 

0030 0002 *m matrix dimensions 

0031 0002 *n matrix dimensions 

0032 0002 *P matrix dimensions 

0033 0000 address of W(K) 

0035 0000 address of V(K) 

0036 0000 address of U(K) 

0037 0000 master counter 
addresses of subroutines 

0046 4000 xstar 

0050 5300 listop 

0052 43080 delay 

0053 4400 mulop 

0054 5200 adop 

0055 4600 matadd 

0056 4700 negmat 

0057 5000 matmul 


0060 to 0077 
.0100 to 0117 


101 


temporary storage 
null matrix 


PROGRAM HYBRIDI CONTINUED 


CELL MACHINE SYMBOLIC REMARKS 

CODE 
0120 0400 *plant transition matrix PHI(1,1) 
0121 0103 PHI(1 ,2) 
0122 0000 PHI(2,1) 
0123 0007 PHI(2,2) 
0140 0236 *plant distribution matrix DEL(1,1) 
0141 0000 
0142 0324 IDED(2 1) 
0143 0000 


0160 to 0177 


0200 


0220 
0221 


0400 


7314 
7650 


filter gain matrix G 
*observable matrix H 


*controller gain matrix AT(1,1) 
ZU) 


0240 to 0257 
0260 to 0277 
0300 to 0317 


noisy observable matrix Z(K) 
filter predicted states XS(K) 
deterministic input matrix 


DI(K-1) 
0320 to 0397 templ 
0340 to 0397 temp2 
0360 to 0377 temp3 
0400 to 0417 temp4 
0420 to 0437 temp5 
0440 to 0457 temp6 
0460 to 0477 temp7 
0500 to 0517 temp8 
0520 to 0537 temp9 
ὑπο 0557 templ0 
0560 to 0577 deterministic input matrix 
DI (K) 
0720 "0252 *filter stable gain matrix G(1,1) 
0721 0000 
0722 0163 G(2,1) 
0723 0000 


4 
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PROGRAM HYBRID I CONTINUED 


CELL 


1000 
1001 
1002 


1003 
1004 
1005 
1006 
1007 
1010 
1011 
1012 
1013 
1014 
1015 
1016 
1017 
1020 
1021 


1022 
1023 
1024 
1025 
1026 


1027 
1030 
1031 
1032 
1033 


1034 
1035 
1036 
1037 
1040 


MACHINE 


CODE 


6205 
2200 
0452 


4025 
6203 
0400 
4025 
2200 
0100 
4205 
2417 
4077 
0400 
4100 
XXXX 
5701 
5477 
6505 


0101 
7057 
0100 
0260 
0260 


0101 
7057 
0100 
0560 
0560 


0101 
7057 
0100 
0300 
0300 


loopl 


C 


SYMBOLIC 
pjf a 
ldc 
0452 
std 25 
ppt ab 
lan 0 
std 25 
ldc 
0100 
sti C 
lcd 17 
std 77 
lan ϐ 
stm 
XXXX 
aob c 
aod 77 
nzb  loopl 
pta 
jpi matmul 
0100 
0260 
0260 
pta 
jpi matmul 
0100 
0560 
0560 
pta 
jpi matmul 
0100 
0300 
0300 
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REMARKS 


for ramp input enter A with 4000 


DA correction for steady state 
error to ramp input 


DA = 0 for step input 


null matrix equal to zero 


-20 


increment address 
increment counter 


initialize XS(K) = 0 


initialize DI(K) = 0 


initialize DI(K-1) = 0 


PROGRAM HYBRID I CONTINUED 


CELL 


1041 
1042 
1043 
1044 
1045 


1046 
1047 


MACHINE 
CODE 


0101 
7055 
0100 
0720 
0160 


0400 
7700 


SYMBOLIC REMARKS 
pta initialize G = Gstable 
jpi matadd 
0100 
0720 
0160 
ldn 0 
hit 
enter A vvith 0000 noise re- 
initialized 
enter A vvith 0001 noise not re- 
initialized 
nzi d 
lice initialization of nose transfers 
0144 noise values from 3400-3543 to 
std 75 3600-3743 so that they can be 
ldc in the same order for each run 
3400 
std 77 
ldc 
3600 
std 76 
ldi 77 
sti 76 
aod 77 
aod 76 
aod 75 
nzb  loop2 
icd 26 set master counter 
std 37 
pta set program loop link 
std 60 
lcd 23 input observables Y(K) by ADC 
std 77 


EE ὃ 


PROGRAM HYBRID I CONTINUED 


CELL MACHINE SYMBOLIC REMARKS 
CODE 

1076 2200 ldc 

1077 1402 1402 analog to digital channel 

1100 4205 stf Ξ 

1101 2200 ldc 

1102 0240 0240 address of Y(K) 

1103 4205 stf f 

1104 7500 loop3 exf select ADC 

1105  xxxxe XXXX 

1106 7600 ina 

1107 4100 stm 

1110  xXExf XXXX 

1111 5704 aob e increment channel, address & 

counter 

1112 2030 ldd 30 

1113 5308 rab f 

1114 5477 aod 77 

1115 6511 nzb  loop3 

1116 2430 led 30 input deterministic signals 

1117 4077 std 77 DI(K) by ADC 

3120 2200 ldc 

(has 1407 1407 analog to digital channel 

1122 4205 stf g 

1:5 -2200 ldc 

1124 0560 0560 address of DI(K) 

1125 4285 stf h 

1126 7500 loop4 ext select ADC 

M27 SRK Gg XXXX 

1130 7600 ina 

1131 4100 stm 

1132 xxx h XXXX 

1133 5704 aob g 

1134 2030 ldd 30 

1185 5908 rab h 

1136 5477 aod 77 

lee 6511 nzb  loop4 


- - ë e ο-- ὃ-- --- -- -ωω ë —— ë ë  ” ” Do  . .. .  .. mo mmo mo umo um 
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PROGRAM HYBRID I CONTINUED 


CELL 


1140 
1141 
1142 
1143 
1144 
1145 


1146 


1147 
1150 
1151 
1152 
1193 
1154 
1155 


1156 
See 
1160 
13851 
1162 
1168 
1164 
1165 
1166 


1167 
1170 


1171 
le 
1173 
1174 
1175 


1176 
ΓΙ; 
1200 


MACHINE 
CODE 


2423 
4062 


“2200 


0240 
4230 
4230 


0101 


7050 
4000 
3600 
3743 
2077 
4100 
3743 


2027 
4204 
0101 
7053 
0022 
XXXX 
6102 
6202 
7700 


2077 
4035 


0101 
7054 
0035 
XXXX 
XXXX 


2030 
9303 
2030 


loops 


SYMBOLIC 
led 23 
std 62 
ldc 
0240 
stf i 
stf j 
pta 
jpi listop 
4000 
3600 
3743 
lidad 7 
stm 
3743 
ldd 27 
stf k 
pta 
jpi mulop 
0022 
XXXX 
nzí 1 
pjf m 
hlt 
ldd 7 7 
std 35 
pta 
jpi adop 
0035 
XXXX 
XXXX 
ldd 30 
rab i 
laa 30° 


REMARKS 


calculate noisy observable Z(K) 


select noise put random number 
in cell 77 


preserve random number list 


calculate V(K) 
scale factor for mulop 


address of sigv 
test mulop overflow 


mulop overflowed, overflow in 
À register 


V(K) = sigv*random number 


Z(K) 7 Y(K) * V(K) 





PROGRAM HYBRID I CONTINUED 
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CELL MACHINE SYMBOLIC REMARKS 
CODE 
1201 5304 rab Τὺ 
1212 5462 aod 62 
]203 Gogo nzb  loop5 
1204 0101 pta calculate W (K) 
1205 7050 jpi listop select noise 
1206 4000 4000 put random number in cell 77 
1207 3600 3600 
1210 3743 3743 
hàll 2077 std 77 preserve random number list 
1212 4100 stm 
1213 3743 3743 
1214 027 Iud 27 scale factor mulop 
1215 4204 stf n 
1216 0101 pta 
1217 7053 jpi mulop 
1220 0021 0021 address of sigw 
1221 XXXX N XXXX 
1222 6102 nzf ο test mulop overflow 
ου 6202 pjf p 
1228 7700 ο hit mulop overflowed, overflow in 
A register 
1] 2257 2077" p ldd 77 
1226 1/033 std 33 W(K) = sigw*random number 
1227 0101 pta DI(K) = -DI(K) 
1230 156 jpi negmat 
1231 0560 0560 
1232 0560 0560 
1233 0101 pta update XS (K) 
1234 7046 jpi xstar 
1235 0101 pta templ = XS(K) - DI(K) 
1236 7055 jpi matadd 
1237 0260 0260 
1240 0560 0560 
1241 0320 0320 


PROGRAM HYBRID I CONTINUED 


CELL MACHINE 
CODE 


1242 
1243 
1244 
1245 
1246 
1247 
1250 
1251 


0101 
7057 
0220 
0320 
0320 
2100 
0320 


SYMBOLIC REMARKS 
pta U(K) = A+| XS(K) - DIK)) 
jpi matmul 

A 

0320 

0320 
ldm 

0320 


4036 address of U(K) 


output U(K) by DAC channel one 


0037 last word output address plus 
one 


update DI(K-1) 
matadd 
0100 
0560 
0300 


internal time delay 
delay for external clock control set 


coarse coarse control equal to 0001 
fine 





PROGRAM HYBRID I CONTINUED 


CELL 


1301 
1302 
1303 
1304 
1305 
1306 


Note! 


MACHINE 
CODE 


9437 
6102 
6202 
7060 
7700 
0036 


un 


SYMBOLIC REMARKS 
aod 37 increment master counter 
net JI 
pjf s 
jpi 60 return to input of Y(K) 
hlt 
0036 first word output address 


If the external clock control is used, the enabling 


interval for ADC sampling must be long enough to allow all 
the samples to be made on each pass. 
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APPENDIX V 
PROGRAM HYBRID II 
V-1 General Description 

Program Hybrid II is a CDC computer program wirtten for the 
control of the simulated hybrid system discussed in Section 16,' 
vshovyn in block diagram form in Figure 16-1.,, and in flow chart 
form in Figure IV-1. 

The program uses a filter to estimate all the states of a 
plant and a controller to synthesize the desired control. The out- 
put of the program is a control for the plant simulated on the PAC 
TR20 analog computer. A stable gain matrix is used for the con- 
troller but a variable gain matrix is computed on line for the 
filter. The program considers all matrices as square matrices and 
will handle up to fourth order systems with minor modification be- 
ginning at cell 2060, and higher order systems with an increase 
in available storage cells for the variables and temporaries. The 
program can be used for plants with only one observable state 
since a matrix inversion process is replaced by division in the 
calculation of the filter gain matrix in subroutine vfgain. The 
following listing of the program is for the second order train plant 
of Figure 4-2 with a sample interval of one second and with a 
pseudo noise to signal power ratio of one taking the noise R(1,1) = 
0.02. An enabling window width of 600 microseconds is sufficient 
to allow all the A/D sampling to be completed on each pass. 

The program uses subroutines, xstar,ligtop, delay, mulop, 
adop, matadd, inegmat, matmul, divop and vfgain, which are 12 
bit variable decimal point routines that are described in Appendix 
VI. 
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For overflow provisions, internal timing control and external 
timing control see Sections 2, 3 and 4 of Appendix III for program 
Optimum Controller. 

V-2 Use of Program Hybrid II 

A. Set up the plant to be simulated on the PACE TR20 analog 
computer similar to Figure 16-1. Connect the three remote control 
connections of the analog computer to those of the ADC so that 
the master clear switch on the CDC 160 digital computer has con- 
trol of the hybrid simulation. 

B. Turn on the ADC/DAC unit, and the CDC 168 arithmetic 
unit. 

C. Load the program at (0000)g and the subroutines at the 
addresses specified by the directory cells. 

D. Enter constants that change with the plant and the order 
of the plant as designated by an asterisk in the following listing 
of the program. For a second order plant enter the initial values 
of P(1,1) in cell 2061 and P(2,2) in cell 2065. For third or higher 
order plants rewrite the section beginning at cell 2060 for the 
initialization of the P matrix. 

E. Set the delay loop or external clock for the sample in- 
terval desired. 

F. Since the gun train plant is a type one system it has a 
steady state error when responding to a ramp input. A DA factor 
was introduced to nullify this steady state error. Hence for ramp 
inputs enter A with (4000) g, run the program from (0000)g and halt 
in cell 1047. Enter A with (0001)g to not reinitialize the noise 
table; otherwise it is restored to the same order, and run the 
program again, In subroutine vfgain to provide for a division 


overflow the variable gain matrix elements are set equal to the 
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stable gain matrix elements and an error flag indicating the 


last countis setin cell 16, 


IL 





PROGRAM HYBRID II 


CELL MACHINE SYMBOLIC REMARKS 
CODE 
0000 7101 jfi 1 run from zero 


0001 2000 2000. location.of program 


constants to be loaded as noted 


0017 0020 counter for null matrix 

0021 0072 *standard deviation of environment 
noise 'sigw' 

0022 0044 *standard deviation of measurement 
noise 'sigv' 

0023 0001 “number of observables 

0024 0004 *produce m*n 

0025 0000 *DA steady state error correction 

0026 0000 number of stages, set equal to 
zero for continuous sampling 

0027 0400 scale factor for mulop, scale 
factor of three 

0030 0002 *m matrix dimensions 

0031 0002 *n matrix dimensions 

0032 0002 *p matrix dimensions 

0033 0000 address of W(K) 

0034 0015 address of sigwsq- variance of 
measurement noise 

0035 0000 address of V(K) 

0036 0000 address of U(K) 

0037 0000 master counter 
addresses of subroutines 

0044 5400 vfgain 

0045 4200 divop 

0046 4000 xstar 

0050 5300 listop 

0052 4300 delay 

0053 4400 mulop 

0054 5200 adop 

0055 4600 matadd 

0056 4700 negmat 

0057 5000 matmul 


PROGRAM HYBRID II CONTINUED 


CELL MACHINE SYMBOLIC REMARKS 
CODE 


0060 to 0077 
0100 to 0117 


0120 
0121 
ο. 
0123 


0140 
0141 
0142 
0143 


0160 to 
0200 
0201 
0202 
0203 


0220 
0221 
0222 
0223 


0240 to 
0260 to 
0300 to 
0320 të 
0340 to 
0360 to 
0400 to 
0420 to 
0440 to 
0460 to 
0500 to 
0520 tọ 
0540 to 
0560 to 


0400 
0103 
0000 
0007 


0236 
0000 
0324 
0000 


ΠΠ 
0400 
0000 
0000 
0000 


7314 
7650 
0000 
0000 


0257 
0277 
0317 
0337 
0357 
0377 
0417 
0437 
0457 
0477 
0517 
0537 
0557 
0577 


temporary storage 
null matrix 


*plant transition matrix PHI(1,1) 
PHI(1,2) 
PHI(2,1) 
PHI(2,2) 
*plant distribution matrix DEL(1,1) 
DEL(2, 1) 


filter gain matrix G 
*observable matrix H 


» "controller gain matrix AT(1,1) 


AT(1,2) 


noisy observable matrix Z(K) 
filter predicted states XS(K) 
deterministic input matrix DI(K-1) 
templ 

temp2 

temp3 

temp4 

temps 

temp6 

temp7 

temp8 

temp9 

templ0 

deterministic input matrix DI (K) 





PROGRAM HYBRID II CONTINUED 


CELL MACHINE SYMBOLIC REMARKS 

CODE 
0600 0005 *environment noise covariance matrix Q(1,1) 
0601 0007 Q(1,2) 
0602 0007 Q(2,1) 
0603 0011 Q(2,2) 


0620 to 0637 


0640 
0641 
0642 
0643 


0660 


0400 
0000 
0103 
0007 


0005 


0700 to 0717 


0720 
0721 
0722 
0723 


0252 
0000 
0163 
0000 


best estimate of error covariance matrix P 


*transpose of plant transition matrix PHIT(1,1) 


PHIT(1,2) 
PHIT(2,1) 
PHIT(2,2) 


*measurement noise covariance matrix R(1,1) 


PNEW 


*filter stable gain matrix G(1,1) 


G(2,1) 


0400 loopl 


4100 


XXXX C 


9701 
5477 
6505 


std 
pjf 
ldn 
std 
ldc 


Str 

lcd 
std 
ldn 
stm 


aob 
aod 
nzb 


0100 
c 

12 
77 

0 


ΧΧΧΧ 
τ 
T 


loopl 


for ramp enter A with 4000 
DA correction for steady state error 


to ramp input 


DA = 0 for step 
set: null matrix equal to zero 


-20 


increment address and counter 
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PROGRAM HYBRID II CONTINUED 


CELL 


2022 
2023 
2024 
2025 
2026 


2027 
2030 
2031 
2032 


2033 


2034 
2035 
2036 
2037 
2040 


2041 
2042 
2043 
2044 
2045 


2046 
2047 
2050 
2051 
2052 


2053 
2054 
2055 
2056 
2057 


SYMBOLIC REMARKS 
pta initialize XS(K) = 0 
jpi matmul 

0100 

0260 

0260 
pta initialize DI(K) = 0 
jpi matmul 

0100 

0560 

0560 
pta initialize DI(K-1) = 0 
jpi matmul 

9100 

0300 

0300 
pta initialize G(K) = 0 
jpi matmul 

0100 

0160 

0160 
pta initialize NEW = 0 
jpi matmul 

0100 

0700 

0700 
pta initialize P(K) = 0 
jpi matmul 

0100 

0620 

0620 
ldc initialize P(0) values must be 


MACHINE 


CODE 
0101 
7057 
0100 
0260 
0260 


0101 
7057 
0100 
0560 
0560 


0101 
7057 
0100 
0300 
0300 


0101 
7057 
0100 
0160 
0160 


0101 
7057 
0100 
0700 
0700 


0101 
7057 
0100 
0620 
0620 


inserted 
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PROGRAM HYBRID Il CONTINUED 


CELL 


2061 
2062 
2063 
2064 
2065 
2066 
2067 


SYMBOLIC REMARKS 
0005 *P(1,1) 
stm 
0620 
ldc 
0400 *P(2,2) 
stm 
0623 Note! this part must be revised 
for third or higher order plants 
pta update filter gain matrix G 
jpi vfgain 
pta update P matrix equal to PNEW 
jpi matadd 
0100 
0700 
0620 
ldn 0 
hit 
enter A with 0000 noise reinitialized 
enter A with 0001 noise not rein- 
itialized 
nef d initialization of noise transfers 
l'ec noise table from 3400-3543 to 
0144 3600-3743 so that they can be in 
std 75 the same order for each run 
ldc 
3400 
std 77 
ldc 
3600 
std 76 
ldi 77 
sti 76 
aod 77 
aod 76 
aod 75 
nzb  loop2 


MACHINE 
CODE 


0005 
4100 
0620 
2200 
0400 
4100 
0623 


loop2 


PROGRAM HYBRID II CONTINUED 


REMARKS 


set master counter 


set program loop link 


input observable Y(K) by ADC 
ADC channel two 


address of Y(K) 


CELL MACHINE SYMBOLIC 
CODE 

2121 2426 d lcd 26 

2122 4037 std 37 

2123 0101 pta 

2124 4060 std 60 

2125 7500 exf 

2126. 1402 1402 

2127 7600 ina 

2130 4100 stm 

2131 0240 0240 

2132 2430 lcd 30 

2138. 4077 std 77 

2134 2200 ldc 

2135 1407 1407 

2136 4205 stf e 

2137 2200 Fac 

2140 0560 0560 

2141 4205 stf f 

2142 7500 loop3 exf 

2143  xxxxe XXXX 

2144 7600 ina 

2145 4100 stm 

2146 xxxx f XXXX 

2147 5704 aob e 

2150. 2030 ldd 30 

7507595304 rab f 

2152 5477 aod 77 

2153 Οὐ] nzb  loop3 

2154 0101 pta 

2155 7050 jpi listop 

2156 4000 4000 

2157 3600 3600 

2160 3743 3743 

2161 2077 ldd 77 
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input deterministic signals DI(K) 
by ADC 


address of DI(K) 


select ADC 


increment channel, address & 
counter 


calculate noisy observable Z(K) 


select noise, put noise in cell 77 


preserve random number list 








PROGRAM HYBRID II CONTINUED 


REMARKS 


calculate V (K) 


scale factor for mulop 


address of sigv 
test mulop overflow 


mulop overflowed, overflow in 
À register 


V(K) = sigv*random number 


Z(K) = Y(K) + V(K) 


La 2 RO OR emm 


CELL MACHINE SYMBOLIC 
CODE 
2162 4100 stm 
28:63 3743 3743 
2164 2027 ldd 27 
2165 4204 stf g 
2166 0101 pta 
2167 7053 jpi mulop 
2170 0022 0022 
Zi] XXXX g XXXX 
2172 6102 ο h 
2173 6002 SE i 
2174 7700 h hlt 
2175 2077 i ldd 77 
2176 4035 std 55 
ZNJ 0101 pta 
2200 7054 jpi adop 
2201 0035 0035 
2202 0240 0240 
2203 0240 0240 
2204 0101 pta 
2.205 7050 jpi listop 
2206 4000 4000 
2207 3600 3600 
2210 3743 3743 
2211 2077 ldd 77 
2212 4100 stm 
211. 3743 3743 
2214 2027 Idd 3 
2215 4204 stf j 
2216 0101 pta 
2217 7053 jpi mulop 
2220 0021 0021 
2221 XXXX j XXXX 
2222 6102 nzi Nk 
2223 6202 pjf ] 


LIS 


calculate VV(K) 
select noise, put noise in cell 77 


preserve random number list 


scale factor for mulop 


address of sigw 


test mulop overflow 


PROGRAM HYBRID II CONTINUED 


CELL MACHINE 


CODE 


SYMBOLIC 


REMARKS 


2224 7700 k hlt 


2225 
2226 


207] 


4033 


mulop overflowed, overflow in 
À register 


33 W(K) = sigw*random number 


templ =[XS(K) - DI(K)] 
matadd 
0260 
0560 
0320 


U(K) = AT [xS(K) - DI w] scalar 
matmul 
0220 
0320 
0320 


0320 
36 





— 


CELL 


2264 
2265 
2266 
2267 


MACHINE 
CODE 


7500 
2411 
7399 
0037 


PROGRAM HYBRID II CONTINUED 


SYMBOLIC REMARKS 
ext output U(K) by DAC channel one 
2411 
out m 
0037 last word output address plus 
one F 
pta update DI(K-1) = DI(K) 
jpi matadd 
0100 
0560 
0300 
pta update filter gain matrix G 
jpi vfgain 
pta update P matrix equal to PNEW 
jpi matadd 
0100 
0700 
0620 
pta internal time delay 
jpi delay for external clock control set 


coarse coarse control equal to 0001, 
fine fine to 0000 


37 increment master counter 

n 

O 

60 return to input of Y(K) 

0036 first word of output address 


RO OO OR emm 


Note! 


If the external clock control is used, the enabling in- 


terval for ADC sampling must be long enough to allow all the 
samples to be made on each pass. 
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APPENDIX VI 

SUBROUTINES FOR CDC 160 DIGITAL PROGRAMS 
VI-1 General Description 

The subroutines described in this appendix handle decimal 
numbers up to twelve bits in length and the format of scaling the 
number is compatible with the CDC 160 computer, the CDC 168 ' 
arithmetic unit and the ADC/DAC conversion unit, 

The following subroutines were used in programs Hybrid I 


and II: 


Memory allocation Subroutine 
3400 to 3543 noise 
4000 to 4107 xstar 
4140 to 4171 scalop 
4200 to 4220 divop 
4300 to 4323 delay 
4400 to 4510 mulop 
4600 to 4654 matadd 
4700 to 4730 negmat 
S000 to 5131 matmul 
5200 to 5252 adop 
5300 to 5370 listop 
5400 to 5531 vígain 


All the subroutines, except mulop and listop which were 
on the library tape, were written to implement programs Optimum 
Controller, Hybrid I and II. Detailed descriptions of their 
purpose, calling sequence, and listing are given on the follow- 
ing pages. 
VI-2 Mulop 

Subroutine mulop uses the CDC 168 arithmetic unit to 
multiply two twelve bit numbers together and postshift the pro- 
duct. The A register is zero on exit for no overflow and is full 
positive or negative on exit depending upon the type of overflow. 
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Overflow can be checked on exit by a non zero jump error halt, 
Mulop uses temporary cells 77 through 74 among which cell 77 
is considered the A register of an arithmetic unit and cell 76, 
the Q register. The execution time is about 400 microseconds 
for no scaling and about 550 microseconds with scaling. 


The calling sequence is: 


std 77 enter multiplier 
pta 
jpi mulop 


address of multiplicand 
scale factor 
ldd 77 load product 


VI-3 Listop 
Subroutine listop shifts a twelve bit word each time it is 
called, lt can be used to shift data ina table in an end off or 
circular mode. Listop uses temporary cells 77 through 75. 
Cell 77 contains the number that is being shifted. 
The calling sequence is: 


pta 
jpi listop 
arg] bit 0: 0 - shift down bit.11:0 - end off 
1 - shift up l-circular 
address of first 
address of last 
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NOISE 


Noise is a list of the first 100 random numbers from CDC 


1604 library subroutine 'RNDEV'. It has a Gaussian distribution 


with essentially zero mean, standard deviation of one and a 


variance of one, Itis stored in cells 3400 through 3543, Pro- 


vision is made in programs so that reinitialization of noise list 


is at discretion of user, 


decimal numbers and the octal numbers scaled by a factor of 


three, 


DECIMAL 


-0.308 
0.963 
05379 
Oe 979 

-0.006 


1.129% 


=e 
SE 
US 
-1.046 
-0.682 
-0.750 
“ 70.421 
^ -0.756 

13369 
ο. 

ZZ 
umo 

0.408 

2.091 
Ue οἱ 
-079S 
-0.531 
-0.156 

1.075 
-0.380 
-0.463 

TUNIS 


OCTAL 


7660 
0366 
0141 
0521 
7776 
7793 
0436 
7304 
7414 
7547 
7364 
7525 
7477 
0153 
7476 
0770 
7729 
1040 
7402 
0150 
1027 
7475 
7414 
7564 
WEE 
0423 
7636 
7611 
0423 


DECIMAL 


-0.654 
ΝΠ 

0.207 
-0.095 
ο... 
-0.312 
-0.725 


>= 003510 


-0 7088 
-0.853 
0.615 
1.155 
0.108 
0.842 
1.442 
1.980 
17367 
2.605 
00; 
ο ο 
0.103 
0.892 
SU 
=0. 960 
allo 
20.25 
ZU DA 
ZI ο 
SI 
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OCTAL 
7530 
7340 
0065 
Lal 
(179 
"7660 
7506 
7646 
7755 
7445 
0239 
0537 
0033 
0327 
0561 
0776 
0536 
1232 
6774 
0040 
0062 
0662 
7979 
7412 
TES 
7477 
7720 
7334 
7476 


DECIMAL 


215315 
20012076 
=O 7160 
-0.507 
-0.840 
2.178 
0.940 
Ue Oo 
0.514 
1:281 
-0.856 
0.174 
29.827 
IO 
0.266 
0.669 
-0.025 
007 2 
aldo 
-2.245 
ο το 
1.284 
1142352 
η 
0.812 
0.543 
-0.491 
-0.081 
0.423 


The following columns are a list of the 


OCTAL 
7656 
7754 
7472 
7575 
7150 
1054 
0360 
0306 
0203 
0510 
7444 
0054 
7454 
pere, 
0104 
0253 
V 
7755 
7303 
6701 
0267 
0510 
0500 
6643 
0317 
0213 
7602 
poe 
0154 





DECIMAL  OCTAL 


ον ου 
-0.463 
d 4075 
-0.476 
0.425 
05019 
0.241 


7636 
7611 
0423 


,7 604 


0154 
0005 
0075 


DECIMAL 


0.043 
-0.386 
-0.814 

0.491 
-1 2998 
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OCTAL 


0012 
7684 
7457 
0175 
7001 


DECIMAL 


E^ 
-0.082 

1.045 
-0.401 


OCTAL 


7033 
7792 
0417 
7631 


SUBROUTINE XSTAR 


Subroutine Xstar solves the controlled recursive filter equation 
XS(K) = (PHI*XS(K-1) + DEL CA! ( XS(K-1)-DI(K-1) #DA JI1-G*H}G*zZ 
and updates the value of XS(K) each time it is called. 

CALLING SEQUENCE 


pta 
jpi Xstar 
CELL MACHINE SYMBOLIC REMARKS 
CODE 
4000 0602 adn 2 return link address 
4001 4061 std 60 
4002 0101 pta 
4003 7057 jpi matmul 
4004 0120 PHI 
4005 260 XS(K-1) 
4006 320 templ 
4007 0101 pta 
4010 7057 jpi matmul 
4011 0160 G 
4012 0200 H 
4013 0340 temp2 
4014 0101 pta 
4015 7056 jpi negmat 
4016 0340 temp2 
4017 0340 -temp2 
4020 0101 pta 
4021 7057 jpi matmul 
4022 0340 temp2 
4023 0320 templ 
4024 0360 temp3 temp3 = -G*H*PHI*XS(K-1) 
4025 0101 pta 
4026 7099 jpi matadd 
4027 0320 templ 
4030 0360 temp3 
4031 0400 temp4 temp4 = [ 1-G*H ]PHI*XS(K-1) 
4032 0101 pta 
4033 7055 jpi matadd 
4034 0260 XS(K-1) 
4035 0300 DI(K-1) 
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| 
| 


SUBROUTINE XSTAR CONTINUED 


CELL 


4036 
4037 
4040 
4041 
4042 
4043 
4044 
4045 
4046 
4047 
4050 
4051 
4052 
4053 
4054 
4055 
4056 
4057 
4060 
4061 
4062 
4063 
4064 
4065 
4066 
4067 
4070 
4071 
4072 
4073 
4074 
4075 
4076 
4077 
4100 
4101 
4102 
4103 
4104 
4105 
4106 
4107 


MACHINE SYMBOLIC 
CODE 


0420 
0101 
7057 
0220 
0420 
0540 
0101 
7054 
0540 
0025 
0540 
0101 
7057 
0140 
0540 
0420 
0101 
7057 
0340 
0420 
0440 
0101 
7055 
0420 
0440 
0460 
0101 
7057 
0160 
0240 
0500 
0101 
7055 
0400 
0460 
0260 
0101 
7055 
0260 
0500 
0260 
7061 


pta 
jpi 


pta 
jpi 


pta 
jpi 


pta 
jpi 


pta 
jpi 


pta 
jpi 


pta 
jpi 


pta 
jpi 


REMARKS 
temp5 temp5 - XS(K-1)-DI(K-1) 


matmul 

AT 

temp5 

templ0 templ0 =A! [ XS(K-1)-DI(K-1) ] 


adop 
templ0 
DA 
templ0 


matmul 

DEL 

templ0 

temp5 temp5 z DELJA IXS(K-1)-DI(K-1) HDA) 


matmul 
temp2 
temps 
temp6 


matadd 

temps 

tempo 

temp7 temp7 = DEL (A! [xs(K-1)-DI(K-1) ] 
+DA ] * [1-G*H] 

matmul 

G 

Z 

temp8 temp8 = G*Z 


matadd 

temp4 

temp7 

XS(K) XS(K)=(PHI*XS (K-1)+DEL[AT[XS (K-1)- 
DI (K- 1) ]+DA [1-G*H] 

matadd 

XS (K) 

temp8 


XS(K) XS(K) = XS(K) + G*Z 


jpi 60 return link 
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SUBROUTINE SCALOP 

Subroutine Scalop shifts the binary point of a number, thereby 
allowing for scaling of octal numbers in the computer on line and 
exits with the result in the A register. In scaling a number down 
the binary point is shifted to the right, however, only left shift 
instructions are available, Hence choose the appropriate number 
of left shifts from the scaling down table. A variable mask de- 
pending upon the type and number of shifts is also given in the 
tables. This mask discards the lower bits. Scalop uses temporary 


cells 77 through 74. 
CALLING SEQUENCE 


pta 
jpi scalop 

number to be scaled 

variable mask 

number of left shifts 
SCALING DOWN TABLE SCALING UP TABLE 
VARIABLE LERT RIGHT l VARIABLE LEFT 
MASKS SHIFTS SHIFTS MASKS SHIFTS 
777 13 l ο Τὸ 1 
0777 12 2 3771 2 
0377 11 3 3770 3 
0177 10 4 3730 4 
0077 7 5 3710 5 
0037 6 6 3700 6 
0017 5 7 3300 7 
0007 4 10 3100 10 
0003 3 11 3000 l1 
0001 D 12 2000 12 
CELL MACHINE SYMBOLIC REMARKS 

CODE 

4140 0602 adn 2 
4141 4077 std Tie 
4142 2177 lai 77 
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SUBROUTINE SCALOP CONTINUED 


CELL 


4143 
4144 
4145 
4146 
4147 
4150 
4151 
4152 
4153 
4154 
4155 
4156 
4157 
4160 
4161 
4162 
4163 
4164 
4165 
4166 
4167 
4170 


4171 


MACHINE 
CODE 


4076 
9477 
2177 
4221 
5427 
2597 7 
4075 
9477 
4074 
2076 
1200 
4000 
4077 
2076 
0102 
4076 
9475 
6504 
2076 
1200 


Xxxx mask 


9077 


7074 


loop 


SYMBOLIC 
std 76 
aod 77 
ldi 77 
stf mask 
aod 77 
lci 77 
std PS 
aod Vel 
std 74 
ldd 76 
lpc 

4000 
std 77 
ldd 76 
ls 1 
std 76 
aod 75 
nzb loop 
ldd 76 
lpc 

XXXX 
rad 77 
jpi 74 


123. 


REMARKS 


number to be scaled 


variable mask 


- number of left shifts 


return link 


sign of number 


increment counter 
shifting completed 


logical product with mask 
scaled number in À register 


at exit 
return main program 


SUBROUTINE DIVOP 


Subroutine Divop forms the quotient of two numbers that are 
equally scaled, 
the product of two numbers, by reversing the arguments for sub- 
routine Mulop, there is no error stop for dividing by zero but 
there are the error stops for full positive and negative overflow. 


Divop uses temporary cells 77 through 73, and directory cell 53 


for mulop. 


CALLING SEQUENCE 


CELL 


4200 
4201 


4202. 


4203 
4204 
4205 
4206 
4207 
4210 
4211 
4212 
4213 
4214 
4215 
4216 
4217 
4220 


Since divop calls subroutine mulop, which forms 





std 77 enter dividend 

pta 

jpi divop 

address of divisor 
scale factor 

ldd 77 quotient 

MACHINE SYMBOLIC REMARKS 
CODE 

0602 adn 2 

4075 std 75 i 
“il 75” ldi 75 address of divisor 

"4074 std 74 

2174 . ldi 74 divisor 

4211 ἅμα: divisor 

5475 aod 75 

2175 ldi 75 scale factor 

4073 std 73 

5475 aod 75 return link 

4206 stf exit 

0101 pta 

7053 jpi mulop 

0073 0073 

XXXx divisor XXXX 

7101 jfi 1 

XXXX exit XXXX 
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SUBROUTINE DELAY 


Subroutine Delay provides a variable-length internal time 


delay loop controlled by the setting of fine and coarse control 


arguments, 


The fine delay loop operates by adding one to the 


complement of the fine control setting and non zero jumping back. 


If the fine control setting is (0000),, then one is added until 


number (7777). negative zero is reached and then the coarse 


control setting is incremented. This count-down of looping 


through the fine control delay and adding one to the complement 


of the coarse control setting is repeated as many times as 


indicated by the coarse delay factor. 


To accurately adjust the internal time delay loops, connect 


a brush recorder to the plant control, and vary the fine and 


coarse delays until the desired sample interval is obtained, 


CALLING SEQUENCE 


pta 
jpi delay 
coarse 
fine 
CELL MACHINE SYMBOLIC REMARKS 
CODE 
4300 0602 adn 2 coarse delay address 
4301 4206 stf a 
4302 0601 adn l fine delay address 
4308 4207 sti. b 
4304 0601 adn 1 return link address 
4305 4214 stf C 
4306 2500 lcm coarse delay complement 
4307  xxxx a XXXX 
4310 4212 stf d 
4311 2500 g lcm ‘fine delay complement 
4312  xxxx b XXXX 
4313 4210 stf ο 
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SUBROUTINE DELAY CONTINUED 


CELL 


4314 
4315 
4316 
4317 
4320 
4321 
4322 
4323 


MACHINE 
CODE 


5607 
6501 
9604 
6506 
7101 
XXXX 
XXXX 
XXXX 


f 


Ω, Ω 


SYMBOLIC 

aof E 

nzb f 

aof d 

nzb g 

jfi l 
XXXX 
XXXX 
XXXX 


REMARKS 
increment fine delay loop 
increment coarse delay loop 
return address 


coarse control 
fine control 





SUBROUTINE MATADD 


Subroutine Matadd takes the sum of two matrices at locations 


A and B, putting the result in location C, and checking for summation 


overflow. Matadd uses temporary cells 77 through 70. 


The di- 


mensions of rows (m) and columns (n) of the matrices are entered 


in cells 30 and 31, respectively, and the address for subroutine 


mulop in directory cell 53. 


CALLING SEQUENCE 


base address of matrix A 
base address of matrix B 
base address of answer matrix C 


mxn 


2 
77 
77 
73 
77 
11 
72 
77 
77 
71 
77 
exit 
30 
77 


mulop 


0031 
0001 
77 
70 
73 


pta 

jpi matadd 

e... MXN 

MATADD = (A) + 
CELL MACHINE SYMBOLIC 

CODE 

4600 0602 adn 
4601 4077 std 
4602 2177 ldi 
4603 4073 std 
4604 5477 aod 
4605 2177 ldi 
4606 4072 std 
4607 5477 aod 
4610 2177 ldi 
4611 4071 std 
4612 9477 aod 
4613 4241 stf 
4614 2030 ldd 
4615 4077 std 
4616 0101 pta 
4617 7053 jpi 
4620 0031 
4621 0001 
4622 2477 lcd 
4623 4070 std 
4624 2173 loop ldi 
4625 3172 adi 
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mxn 
(C) 


REMARKS 


base address of matrix A 


base address of matrix B 


base address of matrix C 


return link 
value of m 
multiplier for mulop 


address of n 
multiply integer, result in 77 


counter = -m*n 
load A(i,j) 
add B(i,j) 


SUBROUTINE MATADD CONTINUED 


CELL 


4626 
4627 
4630 
4631 
4632 
4633 
4634 
4635 
4636 
4637 
4640 
4641 
4642 
4643 
4644 
4645 


4646 


4647 
4650 
4651 
4652 
4653 
4654 


CODE 


4121. 


2173 
6202 
6307 
2172 
6313 
2171 
6211 
0401 
7700 
2172 
6205 
2171 
6303 
0402 
7700 
5473 
5472 
5471 
5470 
6526 
7101 


XXXX exit 


sti 
ldi 
pji 
nji 
ldi 
njf 
ldi 
pji 
ldn 
hit 
ldi 
pji 
ldi 
njf 
ldn 
hlt 
aod 
aod 
aod 
aod 
nzb 
jfi 


MACHINE SYMBOLIC 


71 
73 
a 
b 
We 
ο 
AI 
ο 
] 


2 


Na 


ed 
72 


70 
loop 
1 
XXXX 


REMARKS 


store in C(i,j) 
test for summation overflow 


no overflow, A positive-B negative 
no overflow, A and B positive 


positive overflow, hlt with l in 
A register - 


no overflow, À negative-B positive 


no overflow, A and B negative 


negative overflow, hlt with 2 in 
A register 

increment matrix addresses and 
counter 
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SUBROUTINE NEGMAT 


Subroutine Negmat takes the negative of matrix at location A and 


stores the answer at location B. Negmat uses temporary cells 77 


through 71. 


The dimensions of rows (m) and columns (n) are entered 


into cells 30 and 31 respectively, and the address of subroutine mulop 


in directory cell 53, 


CALLING SEQUENCE 


CELL 


4700 
4701 
4702 
4703 
4704 
4705 
4706 
4707 
4710 
4711 
4712 
4713 
4714 
4715 
4716 
4717 
4720 
4721 
47 22 
4723 
4724 
4725 
4726 
4727 
4730 


pta 
jpi 


negmat 

base address of matrix A 

base address of answer matrix B 
mxn mxn 


NEGMAT of (A)=- (A) 


MACHINE 
CODE 


0602 
4077 
2177 
4073 
5477 
2177 
4072 
5477 
4220 
2030 
4077 
0101 
7053 
0031 
0001 
2477 
4071 
2573 
4172 
5473 
5472 
9471 
6505 
7101 


loop 


XXXX exit 


SYMBOLIC REMARKS 
adn 2 
std 77 
ldi 79 
star 76 base address of matrix A 
aod 77 
ldi 77 | 
std 72 base address of matrix B 
aod 77 
stf exit return link 
ldd 30 value of m 
std 77 multiplier for mulop 
pta 
jpi mulop 
0031 address of n 
0001 multiply integer, result in 77 
σα. vy counter = -m*n 
std 71 
lci 73 load -A(i,j) 
sti 72 store in B(i,j) 
aod 73 increment matrix address and 
aod 72 counter 
aod 71 
nzb loop 
jfi l 
XXXX 
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SUBROUTINE MATMUL 

Subroutine Matmul forms the produce of matrices in locations A 
and B, storing the result in location C, and checking for summation 
and product overflow. Matmul uses temporary cells 77 through 64. 
The scale factor for mulop is entered in cell 27, the dimensions of 
the matrices (m, n, p) in cells 30, 31 and 32, and the address of 
subroutine mulop in directory cell 53, 


CALLING SEQUENCE 


pta 
jpi matmul 
base address of matrix A 
base address of matrix B 
base address of answer matrix C . 
mxp pini m xn 
MATMUL=(A)*(B)=(C) 


CELL MACHINE SYMBOLIC 


REMARKS 
CODE 
5000 0602 adn 2 
5001 4077 std mm 
5002 2177 ldi 77 
5003 4073 std 73 base address of matrix A 
5004 5477 aod 7 
5005. "CAT ldi Th 
5006 4072 std T base address of matrix B 
5007 5477 aod 7 - 
SONJA 72177 ldi 77 
5011 4071 std ΤΙ base address of matrix C 
5012 5477 aod 75 
5013 4244 stf exit 
5014 2027 ldd 27 scale factor for mulop 
5015 4231 stf si 
5016 2032 ldd 32 value of p 
5017 4077 std 77 multiplier for mulop 
5020 0101 pta 
502] 995 jpi mulop - 
5022. 0051 0031 address ofn 








SUBROUTINE MATMUL CONTINUED 


CELL 


9023 
5024 
5025 
5026 
9027 
5030 
5031 
9032 
9033 
5034 
5035 
5036 
9037 
5040 
5041 
5042 
5043 
5044 
5045 
5046 
5047 
5050 
5051 
5052 
5053 
5054 
5055 
5056 
5057 
5060 
5061 
5062 
5063 
5064 
9065 
9066 
5067 
5070 
5071 


MACHINE 
CODE 


0001 
2477 
0601 
4064 
2430 
4065 
2431 
4066 
0400 
4070 
2432 
4067 
2173 
4077 
2072 
4203 
0101 
7053 


loop3 


loop2 


loopl 


xxxx mul 
XXxx sf 


6161 
2077 
4063 
6206 
6314 
0401 
6103 
7101 


XXXX exit 


2070 
6315 
2063 
5070 
6214 
0404 
7700 
2070 
6206 
2063 


a 


b 


SYMBOLIC 
0001 
lcd 77 
adn 1 
std 64 
lcd 30 
std 65 
lcd 3l 
std 66 
ldn 0 
std 70 
lcd oT 
std 67 
ldi 73 
std 77 
ldd 72 
stf mul 
pta 
jpi mulop 
XXXX 
XXXX 
nzf hlt 
ldd 77 
std 63 
pjf a 
njf b 
ldn l 
nzf a 
jfi l 
XXXX 
ldd 70 
njf C 
ldd 63 
rad 70 
pjf d 
ldn 4 
hlt y- m 
ldd 70 
pjf C 
ldd 63 
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REMARKS 


multiply integer, answer in 77 
-p*n 


counter = l-p*n 
-m 

m counter 

-n 

n counter 


set temporary sum equal to Zero 
-p 

p counter 

value of A(i,k) 

multiplier for mulop 

address of B(k,j) 


mulop overflowed 
a(i,k)*B(k,j) 


check summation overflow 
i.e. A(1,1)*B(1,1)+ 
A(1,2)*B(2,1) 
absolute jump forward 


return address 

temporary sum 

no overflow possible, opposite 
signs 

check temporary sum 


positive summation overflow, hlt 


...With 4 in A register 


. no overflow possible, opposite 


signs 


SUBROUTINE MATMUL CONTINUED 





CELL MACHINE SYMBOLIC REMARKS 

CODE 
5072" 5070 rad 70 check temporary sum 
5073 6305 njf d 
5074 0405 ldn 5 negative summation overflow, hlt 
5075 7700 hlt with 5 in A register 
5076 2063c ldd 63 
5077 5070 rad 70 temporary sum 
5100  2031d ldd 3n increment B(k,j) by n 
5101 5072 rad 72 
5102 0401 ldn 1 increment A(i,j) by one 
5103 -5073 rad 73 | 
5104 5467 aod 67 loop for p products 
5105 6546 nzb loopl 
5106 2070 ldd 70 temporary sum 
5107 4171 sti AI store in C(i,j) 
5110 0401 ldn aud 
5111 5071 rad i increment C(i,j) 
5112 2432 led 32 -p 
5113 5073 rad 73 return A(i,k) to base address 
5114 2064 ldd 64 
5115 5072 rad 72 return B(k,j) to base address 
5116 5466 aod 66 plus one loop for n columns 
5117 254 nzb laop2 Eo 
5120 2032 ldd 32 p 
5121 5073 rad 73 advance address of A(i,k) to 

a im 3 A(i,k) τα 

5122 124818 a icd 31l -n 
5123 50%2 rad 72 return B(k,j) to base address 
5124 5465 aod 65 loop of m rows 
5125 6574 nzb loop3 | 
5126 6650 pjb exit absolute jump backward to exit | 
odo. CD njb exit | 
5130 0403 hlt ldn 3 hlt vvith 3 in A register, mulop 
5131 7700 hlt overflowed, check scale factor 
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cell 27 


SUBROUTINE ADOP 


Subroutine Adop takes the sum of scalars at addresses A and 


B and puts the result in C, and checks for summation overflow. 


The A register is zero on exit for no overflow and is full positive 


or negative on exit depending upon the type of overflow. Over- 


flow can be checked on exit by a non zero jump error halt. Adop 


uses cells 77 through 73. 


CALLING SEQUENCE 


CELL 


5200 
5201 
5202 
9203 
5204 
9205 
9206 
5207 
9210 
5211 
5212 
5213 
5214 
9215 
9216 
5 215 
9220 
9221 
92422 
5223 
5224 
9225 
5226 


pta 
jpi adop 

address of A 

address of B 

address of C 
MACHINE SYMBOLIC 
CODE 
0602 adn 2 
4077 std 77 
2177 ldi 77 
4076 std 76 
5477 aod 77 
2177 ldi 71 
4075 std 19 
5477 aod 77 
2177 ldi QI 
4074 std 74 
5477 aod 77 
4073 std 73 
2176 ldi 76 
6202 pjf a 
6314 njf b 
2175 ldi 75 
6327 njf C 
3176 adi 76 
4174 sti 74 
6226 pjf d 
6303 njf errorP 
2073 ldd 73 
4020 std 20 
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REMARKS 


address of A 


address of B 


address of C 


return link 
value of A 


value of B 

no overflow possible, opposite 
signs 

C=A+B 

no overtlow 


SUBROUTINE ADOP CONTINUED 





CELL MACHINE SYMBOLIC REMARKS 
CODE 
522Z 2216 errorP ldi fullP 
5230 4077 std 77 
σοι 6221 pjf exit positive summation overflow 
DZ tege, | πο Ὦ ldi 7:5 
9233 6214 pjf c 
5234 3176 ' adi 76 
5235 4174 sti 74 C=A+B 
Se EN njf d no overflow possible, opposite 
5237 5203 pji error N signs 
5240 2073 ldd το 
5241 4020 std 20 
5242 2204 errorN ldi fullN 
5243 4077 std 77 
5244 6306 njf exit negative summation overflow 
5245 3777 fullP 
5246 4000 fullN 
5247 3176c adi 76 
5250 4174 sti 74 C=A+B | 
9251 0400 d ldn 0 
5252 7073 exit  jpi 73 
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SUBROUTINE VFGAIN 

Subroutine Vfgain solves the following recursive relationships 
in updating the filter gain matrix G, and the conditional covariance 
matrix P. 

G(K) = P(K)*HT[H*P(K)*HT + RJ 71 

P(K) = PHI*[P(K-1) - G(K-1)*H*P(K-1)]*PHIT + Q 
CALLING SEQUENCE 


pta 
jpi vigain 
CELL MACHINE SYMBOLIC REMARKS 
CODE 
5400 0602 adn 2 
5401 4061 std 60 return link 
5402 0101 pta 
5403 7057 jpi matmul 
5404 0620 P 
5405 0200 ut 
5406 0320 templ templ = P*HÎ 
5407 0101 pta 
5410 7057 jpi matmul 
5411 0200 H 
5412 0320 templ 
5413 0340 temp2 
5414 0101 pta 
5415 7055 jpi matadd 
5416 0340 temp2 
5417 0660 R 
5420 0360 temp3 temp3 = H*P*HÎ + R scalar 
5421 2430 lcd 30 -m 
5422 4063 std 63 
5423 2200 ldc since there is only one observ- 
5424 0160 0160 able for this plant, matrix 
5425 4233 stf a inversion can be replaced by 
5426 2200 ldc division to obtain the filter 
5427 0320 0320 gain matrix elements 
5430 4205 stf b 
5431 2200 ldc 
5432 0360 0360 
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SUBROUTINE VFGAIN CONTINUED 


CELL MACHINE SYMBOLIC REMARKS 
CODE 
5433 4210 stf © 
5434 2100 loopl ldm 
5435 ΧΧΧΧῸ XXXX 
5436 4077 std 77 dividend for divop 
9437 2027 ldd Qs scale factor for divop 
5440 4204 stf sf 
5441 0101 pta 
5442 7045 jpi divop 
5443 XXXXC xxxx  _ G(1,1)=P(1,1)*H!/(H*P*H!+R) 
5444 xxxx sf xxxx  G(2,1)=P(2,1)*H!/(H*P*H! +R) 
5445 6011 zjf d 
5446 2037 ldd 9 master counter 
5447 4016 std 16 set error flag for divop overflow 
5450 0101 pta 
5451 7055 jpi matadd set G=G stable for divop overflow 
5452 0100 zero 
5453 0720 Gstable 
5454 0160 G 
5455 7061 jpi 61 return main program 
5456 2077 ldd 77 quotient from divop 
5457 4100 stm 
5460 xxxxa XXXX store G matrix elements 
5461 2030 ldd 30 
5462 5302 rab 2 increment addresses of G and 
templ 
5463 2030 ldd 30 
5464 5327 rab b 
5465 5463 aod 63 increment counter 
5467 0101 pta 
5470 7057 jpi matmul update P matrix 
5471 0200 H 
5472 0620 P 
5473 0400 temp4 
5474 0101 pta 
5475 7057 jpi matmul 
5476 0160 G 
5477 0400 temp4 
5500 0420 temp5 
5501 0101 pta 
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SUBROUTINE VFGAIN CONTINUED 


CELL 


9902 
9503 
9504 
5905 
5906 
9907 
5510 
5511 
9912 
5913 
9514 
9515 
5916 
9917 
5520 
9521 
9922 
9923 
5924 
5525 
9926 
5927 
5930 
9931 


MACHINE SYMBOLIC 


CODE 


7056 
0420 
0420 
0101 
7055 
0620 
0420 
0440 
0101 
7057 
0120 
0440 
0460 
0101 
7057 
0460 
0640 
0700 
0101 
7055 
0700 
0600 
0700 
7061 


ipi 
pta 


jpi 


pta 
jpi 


pta 
jpi 


pta 
jpi 


jpi 


negmat 
tempo 
-tempo 


matadd 
E 
tempo 
temp6 


matmul 
PHI 
temp6 
temp7 


matmul 
temp7 
PHIT 
PNEW 


matadd 
PNEW 
Q 
PNEW 
61 
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REMARKS 


tempS = - G*H*P 


PNEW=PHI(P - G*H*P)PHI? 


PNEW = PNEW + Q 
return main program 
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